# Mechanical strain induced topological phase changes of few layer ZrTe5

Zoltán Tajkov  HUN-REN Centre for Energy Research, Institute of Technical
Physics and Materials Science, 1121 Budapest, Hungary  Centre of Low
Temperature Physics, Institute of Experimental Physics, Slovak Academy of
Sciences, Košice SK-04001, Slovakia<EMAIL_ADDRESS>Konrád Kandrai
HUN-REN Centre for Energy Research, Institute of Technical Physics and
Materials Science, 1121 Budapest, Hungary  Dániel Nagy  ELTE Eötvös Loránd
University, Department of Physics of Complex Systems, 1117 Budapest, Hungary
Levente Tapasztó  HUN-REN Centre for Energy Research, Institute of Technical
Physics and Materials Science, 1121 Budapest, Hungary  János Koltai  ELTE
Eötvös Loránd University, Department of Biological Physics, 1117 Budapest,
Hungary  Péter Nemes-Incze  HUN-REN Centre for Energy Research, Institute of
Technical Physics and Materials Science, 1121 Budapest, Hungary
<EMAIL_ADDRESS>
###### Abstract

Understanding the topological aspects of the band structure of solids has
fundamentally changed our appreciation of their properties. The layered, van
der Waals transition-metal pentatelluride ZrTe5 has proven on numerous
occasions to be an excellent candidate for the study of controllable
topological phase transitions. Here, we investigate the topological phase
diagrams of monolayer and bilayer forms of ZrTe5, under mechanical
deformations using ab initio techniques. We find that mechanical deformation
can close the monolayer’s topological gap, while the bilayer exhibits richer
phase diagram, including both topological insulating, trivial metallic and
insulating phases. The bilayer is predicted to be on the topological phase
boundary. We also address the preparation of monolayers, using _ab initio_
simulations and experimental scanning tunneling microscopy measurements. We
confirm that while monolayer ZrTe5 is difficult to exfoliate without
compromising its crystalline structure, bilayers offer a more stable
alternative, revealing the complexities and limitations of using gold
substrates for monolayer exfoliation.

topological insulators, monolayer, DFT, phase diagram, SIESTA, VASP

††preprint: APS/123-QED

##  I Introduction

Since the seminal discovery of time-reversal invariant topological insulators,
a controllable topological phase transition in two-dimensional materials
remains a targeted objective in the field of solid-state physics [1, 2, 3, 4].
Such phase transitions are characterized by changes in the topological
properties of the materials band structure [5, 6]. One effective approach for
accomplishing this involves tuning the lattice parameters of topological
materials, potentially through mechanical deformation. The transition-metal
pentatelluride, ZrTe5, resides at the boundary between weak and strong
topological phases [7, 8, 9, 10, 11, 4, 12, 13], thus enabling the use of
strain to toggle between these phases. Indeed, transitions in bulk or thick
exfoliated crystal samples have been observed in prior work [14, 15, 16, 3, 4,
13]. The material also exhibits a variety of other noteworthy phenomena. These
include the planar Hall effect [17], anomalous Hall effect [18], and chiral
magnetic effect [19]. For bulk samples, a Dirac semimetal-like band structure
has been predicted based on early studies of quantum oscillations [20],
optical conductivity [21], and photoemission [19].

Significant effort has been devoted to and still focuses on understanding the
properties of bulk ZrTe5. Yet, studies on few-layer to monolayer crystals are
not as comprehensive, especially concerning the behavior of the ZrTe5 system
under mechanical deformations. The bulk crystal is known to be highly
sensitive to minor changes in its lattice parameters. Transport measurements
reveal non-monotonic strain dependence and a switch in the sign of the Dirac
mass with deformation [9]. Existing ARPES measurements show that under
external strain, the band gap increases with tensile strain and decreases with
compressive strain [22].

Intentionally induced, controllable mechanical distortion effectively switches
topological phases in the crystal. This can significantly change the
electrical properties of a material, including its topological characteristics
[23, 24, 25, 26, 27, 13]. In contrast, for thin [28, 29] to monolayer [30]
samples, experimental studies of ZrTe5 or its sister compound HfTe5 remain
scarce. Nonetheless, recent advancements have been made in the production of
thin ZrTe5 samples through bottom-up growth [31] and Al2O3 assisted
exfoliation [30].

Theoretical studies indicate that the monolayer is a $\mathbb{Z}_{2}$ time-
reversal invariant topological insulator [32]. In this work, we investigate
the topological phase diagram of monolayer and bilayer ZrTe5 crystals under
mechanical distortions, employing ab initio techniques. We also aim to
experimentally address the gap in the preparation of mono and few-layer
samples by exploring the exfoliation of thin ZrTe5 layers onto gold
substrates, following established methods [13]. Exfoliation on gold is a well
known technique for preparing monolayers of van der Waals (vdW) materials that
possess chalcogen terminations on their surfaces [33]. Our results reveal that
this method is unviable for ZrTe5, as the strong attraction between the first
ZrTe5 layer and the gold substrate compromises the crystalline structure of
the initial layer.

##  II Results and Discussion

###  II.1 Topological phase diagrams of monolayer and bilayer ZrTe5

Zirconium pentatelluride (ZrTe5) adopts a layered orthorhombic crystal
structure. Layers spanning the $x-y$ plane form trigonal prismatic ZrTe3
lines, which are connected by zig-zag chains of Te atoms, as illustrated in
Fig. 1a). The chains form two-dimensional sheets that are stacked along the
$z$-axis, creating a layered structure. The primitive cell of bulk ZrTe5
contains two Zr and ten Te atoms [32]. Each ZrTe5 layer is nominally charge
neutral, and the interlayer coupling is of vdW type [32]. The mono- and
bilayer initial geometry was created by simply cutting one or two layers from
the bulk. The initial geometry underwent full relaxation via ab initio
calculations (see Methods section for further details). The starting point for
the relaxation were the lattice constants derived from our previous X-ray
diffraction measurements of the bulk sample [13]. The relaxed and experimental
geometry parameters are listed in Table 1. Our analysis reveals that the in-
plane unit cell parameters and internal atomic coordinates closely align with
their 3D bulk experimental values, corroborating previous findings [32]. The
main difference between the density functional theory (DFT) codes lies in
their prediction of interlayer distances in the bilayer. In particular, the
SIESTA geometry relaxations yield a distance that is 0.13 Å less than that
predicted by VASP. Further discussion on this discrepancy will be provided
later.

|  |  $a$ [Å] |  $b$ [Å] |  $c$ [Å]
---|---|---|---|---
Experimental | Bulk | 1.994 | 7.265 | 13.724
SIESTA | Monolayer | 1.999 | N.A. | 13.870
Bilayer | 1.999 | 7.172 | 13.886
VASP | Monolayer | 2.002 | N.A. | 13.773
Bilayer | 2.001 | 7.307 | 13.758
Table 1: Unit cell dimensions for ZrTe5 as derived from the analysis of X-ray
diffraction data at 300 K [13] and from relaxed DFT calculations. We note that
in the few layer cases the $b$ unit cell vector represents the interlayer
distance.

It has been previously predicted that a single layer of ZrTe5 is a
$\mathbb{Z}_{2}$ time-reversal invariant topological insulator with a band gap
of approximately $100\,\mathrm{meV}$ [32]. In Fig. 1 b)-c), we present the
band structures of the unperturbed crystals, computed by both the SIESTA and
VASP implementations of DFT. The blue solid lines represent the monolayer
case, while the red dashed lines correspond to the bilayer. In the case of the
monolayer, both SIESTA and VASP predict that the undistorted crystal is a
topological insulator, albeit with slightly different band gaps: 110 meV for
VASP and 70 meV for SIESTA. This consensus is broken in the bilayer scenario.
VASP predicts that the bilayer is a trivial insulator with a band gap of 70
meV, whereas SIESTA anticipates a metallic phase.

Figure 1:  a) The crystal structure of the bilayer. b) The band structures of
the undistorted crystal obtained by SIESTA. The blue solid lines (dashed red)
correspond to the monolayer (bilayer) case. c) The band structures of the
undistorted crystal obtained by VASP. The blue solid lines (dashed red)
correspond to the monolayer (bilayer) case. d)-g) The phase diagrams of the
electronic structure of the crystal under mechanical strain obtained by two
different DFT codes for the mono- and the bilayer cases. At every point the
size of the gap was calculated and a sign has been assigned to it according to
the topological favor of the gap. The negative gap corresponds to the
topological insulating phase, while the positive gap to the trivial phase.

As previously demonstrated, the bulk ZrTe5 is highly sensitive to minor
deformations in the crystal structure [13]. This sensitivity is also evident
in the case of few-layer crystals. In Fig. 1 d)-g), we present the topological
phase diagrams for mono- and bilayer ZrTe5 under in-plane mechanical strain,
as a function of the strain along the $x$ and $y$ directions:
$\varepsilon_{xx}$ and $\varepsilon_{yy}$. For each point on the diagram, the
size of the gap is calculated, and a sign is assigned to it to indicate its
topological character. A negative gap corresponds to the topological
insulating phase, whereas a positive gap signifies the trivial phase. Our
exploration involved analyzing the crystals using both SIESTA and VASP. As
illustrated in Fig. 1 d) and f), SIESTA and VASP concur that the monolayer
resides at the boundary between topological and trivial (semimetallic) phases.
Applying a modest uniaxial strain, the system can be tuned from a topological
insulating to metallic phase. These data suggest that the monolayer inherits
from the bulk crystal the property of residing at the boundary of a
topological phase transition.

In Fig. 1 e) and f), we compare the topological phase diagrams of bilayer
ZrTe5. SIESTA predicts that the bilayer resides at the boundary between the
topological insulator and metallic phases. The diagrams indicate that applying
a moderate biaxial strain can transition the system into a topological phase,
while an experimentally achievable biaxial strain can induce a trivial
insulating phase through the metallic regime.

The differences seen in the band structures of the bilayer are also evident
when comparing the phase diagrams. In Fig. 1 f), it can be observed that,
according to the VASP calculations, the bilayer resides further away from the
phase boundary and is in a topologically trivial phase with a significant band
gap of 70 meV. However, the central feature from SIESTA, that biaxial strain
can be used to transition the system into the topological phase, is confirmed
by VASP, but much larger strain is required to realize the transition.

Figure 2:  The topological phase diagram of the bilayer ZrTe5 calculated by
VASP. The interlayer distance has been matched to the SIESTA optimized
interlayer distance.

Further investigation into the role of geometry and lattice distortions
clarifies the origin of the differences between VASP and SIESTA in the bilayer
case. In Fig. 2, we present the topological phase diagrams of the bilayer
calculated with VASP. This time, however, the interlayer distance has been
modified to align with the SIESTA optimization results; specifically, we
adjusted the distance between the layers by making it 0.13 Å smaller. With
this adjustment the predictions of the two codes converge substantially.
Notably, our choice to align VASP with SIESTA in this aspect was driven by our
prior observations regarding the accuracy of SIESTA’s predictions. These
observations, as discussed in Tajkov et al. (2022), revealed that SIESTA
consistently yielded results that demonstrated better agreements with
experimental data, particularly in the context of bulk properties of ZrTe5
[13].

###  II.2 Exfoliation of thin ZrTe5 onto a gold substrate

Preparing large monolayers of chalcogen terminated van der Waals materials can
be easily achieved by exfoliating them onto clean gold surfaces [33]. This has
been tested a wide range of dichalcogenides, including MoS2, WSe2, WS2, etc.
and importantly for Te terminated crystals [33, 34]. For ZrTe5, an Al2O3-based
exfoliation method has been employed [30], yielding large aspect ratio
monolayer flakes that have a lateral size of $\sim$10 $\upmu\mathrm{m}$.

To prepare the sought after large gap 2D topological insulator phase of ZrTe5,
we have explored the exfoliation of ZrTe5 onto gold substrates. This method
yields thick crystals, which show a good crystalline order on the top layer,
as measured by scanning tunneling microscopy (STM) in our previous work [13].
The underlying principle for gold-assisted exfoliation of chalcogen-terminated
materials is that the layer in direct contact with the gold forms a stronger
adhesion to the gold surface than it does to the bulk crystal itself[33, 34].
This approach results in cm-sized crystals in the case of MoS2 [34].

In our exfoliation attempts, we employ the well-known technique of evaporating
gold onto a mica surface, which predominantly yields $[111]$ facets. The gold
is subsequently stripped from the mica [33], and the clean, mica-facing side
is used for the exfoliation of ZrTe5. An optical microscopy image of such an
exfoliation attempt is presented in Fig. 3a. In these experiments, the
faintest optical contrast corresponds to the monolayer flakes. In our
experiments, these areas are indicated by red arrows. All exfoliation
experiments are conducted in a glove-box environment to minimize the risk of
ZrTe5 oxidation. We further investigate the faint monolayer areas using low-
temperature (9K) STM measurements. An exemplar of these measurements is
displayed in Fig. 3b. The results reveal a highly disordered structure, making
it difficult to distinguish the characteristic Zr lines of ZrTe5. This
disorder is further confirmed in the Fourier transform of the image, as shown
in Fig. 3c. In contrast, STM measurements on the top layer of a multi-layer
flake yield atomic resolution, and the Fourier transform of the image reveals
the periodicity associated with the Zr chains, as evident in Fig. 3d, e.

We have consistently observed the disordered nature of the monolayer candidate
flakes across many tens of samples; none of these display clear atomic
resolution in STM images. Given these challenges, it may be beneficial to
consider alternative fabrication methods, such as the growth of monolayer
ZrTe5 via molecular beam epitaxy. In contrast, thicker flakes invariably
permit atomic resolution. To shed light on this puzzling observation, we
turned to _ab initio_ calculations to model the interaction between ZrTe5 and
the gold $[111]$ surface. We performed the computational study of the
interface between the gold surface and multi-layer ZrTe5 using the VASP code
exclusively, as its capabilities are well-suited for such calculations.
Further details can be found in the Method section. Example results of the
relaxed geometry of the monolayer and the bilayer on top of a gold slab can be
seen in Fig. 3f, g. These show that the monolayer becomes highly distorted on
top of the gold surface.

Figure 3:  a Optical microscopy image of exfoliated ZrTe5 crystals on a gold
surface. The thinnest crystals on the surface are marked by red arrows. b STM
topographic image in the thinnest area, as marked by red arrows in a). The
image shows a disordered structure. c Fourier transform image of the data in
b). d STM topographic image of a few layer ZrTe5 flake, showing high
crystallinity, as shown by the Fourier transform in e). Inset: atomic
resolution STM image of the surface of a thick ZrTe5 flake, red bar marks the
unit cell period in the _c_ direction. e Fourier transform of the ZrTe5 area
of d). Red dashed circles mark the Fourier component of the periodicity in the
_c_ direction. f-g The relaxed geometry obtained by the VASP calculation in
the monolayer and bilayer case.  Figure 4:  Distribution of atomic distances
in the few layer ZrTe5 crystals on gold [111] surface in different
configurations as calculated by the VASP code. a-c: Comparing the atomic
distances in the ZrTe5 monolayer for the freestanding case (blue) and relaxed
on top of gold $[111]$ (orange). d-f For bilayer ZrTe5, we compare atomic
distances between the ”top” layer, which is not in contact with the gold
surface, and the ”bottom” layer, which is in direct contact with the gold
surface. g Comparing the Te-Te distances of the ”top” layer on gold to the
freestanding bilayer. f The Te-Te distances of a bottom layer of a bilayer on
gold and that of a monolayer. i Showing the distortion of the Au surface
itself.

To quantify the distortion of the ZrTe5 layer when in contact with the gold
surface, we provide histograms of the atomic distances. In Fig. 4, we display
the interatomic distances in ZrTe5 mono- and bilayers, both in gold-surface-
adhered and freestanding configurations. The first row focuses specifically on
the monolayer case. Panels a)-c) reveal that, in this monolayer configuration,
interatomic distances undergo notable alterations, leading to a significantly
deformed crystal structure compared to the freestanding case. This change is
particularly prominent in the Te-Te distances, as highlighted in panel c).

Transitioning to the bilayer case, observations indicate that in this two-
layer arrangement, the layer proximal to the gold surface (labelled ”bottom”
on Fig. 4) exhibits characteristics similar to a monolayer. Meanwhile, the top
layer, aside from minor variations, largely remains intact. The most striking
evidence of this phenomenon is apparent in the alterations to the Te-Te atomic
distances, as seen in panel f. The relative stability of the top layer is
apparent when comparing its atomic distances to those of a freestanding
monolayer, as illustrated in panel g). Moreover, panel i) reveals that the
interatomic distances within the gold surface also vary depending on the
presence or absence of ZrTe5 on the substrate.

Based on these findings, we conclude that exfoliating a single layer of ZrTe5
onto a gold surface induces substantial disorder in its crystalline order. In
contrast, for bilayer configurations, the layer in contact with the gold
undergoes distortion while the top layer largely remains intact. It is
important to mention that future research should explore how closely the
electrical properties of the top layer in this distorted bilayer resemble
those of a freestanding monolayer.

To gain a more comprehensive understanding of the interaction between the gold
surface and few-layer ZrTe5, we examined the bonding energies of the
heterostructure. The bonding energy between the gold slab and a monolayer of
ZrTe5 is 2.43 eV per unit cell, which is more than double the bonding energy
between a five-layer ZrTe5 stack and a monolayer, recorded at 1.10 eV per unit
cell. Additionally, we calculated the energy required to separate the top
layer from a bilayer exfoliated onto a gold surface, finding it to be 1.16 eV
per unit cell. See section IV for more details. These results indicate that
while the exfoliation of monolayers is highly likely, their crystalline
structure becomes compromised. This is corroborated by our experimental
observations, which reveal expansive areas in optical microscopy with the
weakest optical contrast, suggesting the presence of monolayers. However, STM
measurements indicate that these areas are disordered. Moreover, bilayers are
marginally more likely to form than monolayers through the exfoliation
process, due to the 60 meV stronger bonding discussed above. The validation of
this prediction through experimental methods is a subject of ongoing
investigation.

##  III Conclusions

We have probed the topological phase diagram of both single- and double-layer
standalone ZrTe5 under mechanical stress, utilizing ab initio methods. Our
investigation unveiled that the monolayer exhibits a topological insulating
phase. Under the influence of modest uni- or biaxial strain, this phase can
transform into a metallic state, offering valuable flexibility in tuning its
properties. Unlike the monolayer, the bilayer crystal has areas in its phase
diagram where it becomes a trivial insulator. Moreover, it can be tuned into
either metallic or trivial insulating phases as well. Upon closer examination,
it becomes evident that when a monolayer of ZrTe5 is exfoliated onto a gold
surface, it exhibits a distinctive structure, adding complexity to the
exfoliation process. This observation is supported by our experimental
results.

##  IV Methods

The optimized geometry and electronic properties of the crystal were obtained
by the SIESTA and VASP implementations of density functional theory (DFT) [35,
36, 37, 38, 39, 40].

SIESTA

SIESTA employs norm-conserving pseudopotentials to account for the core
electrons and linear combination of atomic orbitals to construct the valence
states. The generalised gradient approximation of the exchange and the
correlation functional was used with Perdew–Burke–Ernzerhof parametrisation
[41] and the pseudopotentials optimised by Rivero et al. [42] with a
double-$\zeta$ polarised basis set and a realspace grid defined with an
equivalent energy cutoff of 350 Ry for the relaxation phase and 900 Ry for the
single-point calculations. The Brillouin zone integration was sampled by a
30$\times$18$\times$1 Monkhorst–Pack $k$-grid for both the relaxation and the
single-point calculations [43]. The geometry optimisations were performed
until the forces were smaller than 0.1 eV nm-1 . The choice of
pseudopotentials optimised by Rivero et al. ensures that both the obtained
geometrical structures and the electronic band properties are reliable. After
the successful self consistent cycles the necessary information was obtained
by the sisl tool [44]. The spin orbit coupling was taken into account in the
single point calculations. In the SIESTA case we have calculated the
$\mathbb{Z}_{2}$ invariant using a home made tool based on the numerical
method developed by Fukui and Hatsugai [5, 45].

VASP

The VASP code was used with projector-augmented wave pseudopotentials and the
optB86b-vdW functional throughout our computation which has a non-local
correlation correction that approximately accounts for dispersion interactions
[46]. We have found that employing this functional gives the most accurate
geometries in comparison with experiments. The plane-wave cutoff energy was
set to $500\,\rm eV$ in all calculations. The structural optimizations were
performed with $20\times 6\times 1$ Monkhorst-Pack set of $k$-points until all
atomic forces fell below $3\,\rm meV/$Å. The band structure calculations were
performed in the usual way, first the charge density was obtained in a self-
consistent non-collinear calculation with a 24$\times$8$\times$1 Monkhorst-
Pack set, and in the second step the band structure was evaluated along the
high-symmetry lines of the Brillouin zone.

VASP calculations of ZrTe5 on gold

We simulated the gold (111) surface by taking 6 layers of gold and fixed the
bottom 2 layers during optimization. In order to get a commensurate super cell
with the ZrTe5 we needed 7 unit cells of the crystal. The supercell unit
vectors are the following: $\mathbf{A}_{1}=-3\mathbf{a}+\mathbf{c}$ and
$\mathbf{A}_{1}=4\mathbf{a}+\mathbf{c}$, where $\mathbf{a}=\left(a,0,0\right)$
and $\mathbf{c}=\left(0,0,c\right)$. This gives a 27.192° rotation between the
gold and the ZrTe5 unit cell vectors and a slight, $<2\%$ mismatch between the
ZrTe5 and the gold slab. Altogether the gold slab + monolayer ZrTe5 supercell
contains 312 gold atoms, 14 Zr atoms and 70 Te atoms. These calculations were
performed using a single Gamma-point, which is sufficient for large unit cells
like this.

STM measurements

The ZrTe5 crystals, sourced from hqgraphene.com, were exfoliated onto gold
substrates within an inert glovebox setting. These samples were then moved to
the UHV STM chamber, through a vacuum shuttle. We conducted STM measurements
using an RHK PanScan Freedom microscope, operating at both 300 K and 9 K in an
ultra-high vacuum with a base pressure of 5$\times$10-11 Torr. For the STM
tips, mechanically cut Pt/Ir (90%/10%) wire was employed. A high-magnification
optical microscope facilitated the precise positioning of the STM near the
chosen ZrTe5 crystals.

## Competing interests

The authors declare no competing interests.

## Data Availability

Any further data are available from the corresponding author upon reasonable
request.

###### Acknowledgements.

The work was conducted within the framework of the Topology in Nanomaterials
Lendulet project, Grant No. LP2017-9/2017, with support from the European
H2020 GrapheneCore3 Project No. 881603 and the Korea-Hungary bilateral project
(TT-IPARI-KR-2022-00006). Financial support from Élvonal Grant KKP 138144,
NKFIH OTKA grant K132869, K134437 and TKP2021-NKTA-05 grant is also
acknowledged. PNI acknowledges the support of the Janos Bolyai Research
Scholarship of the Hungarian Academy of Sciences. JK acknowledges that this
research was supported by the Ministry of Culture and Innovation and the
National Research, Development and Innovation Office within the Quantum
Information National Laboratory of Hungary (Grant No.
2022-2.1.1-NL-2022-00004). ZT acknowledges financial support from Slovak
Academy of Sciences project IMPULZ IM-2021-42.

## References

  * Murakami  _et al._ [2007] S. Murakami, S. Iso, Y. Avishai, M. Onoda, and N. Nagaosa, Tuning phase transition between quantum spin Hall and ordinary insulating phases, Physical Review B 76, 205304 (2007).
  * Murakami [2007] S. Murakami, Phase transition between the quantum spin Hall and insulator phases in 3D: emergence of a topological gapless phase, New Journal of Physics 9, 356 (2007).
  * Lu  _et al._ [2017] J. Lu, G. Zheng, X. Zhu, W. Ning, H. Zhang, J. Yang, H. Du, K. Yang, H. Lu, Y. Zhang, _et al._ , Thickness-tuned transition of band topology in ZrTe5 nanosheets, Physical Review B 95, 125135 (2017).
  * Mutch  _et al._ [2019] J. Mutch, W.-C. Chen, P. Went, T. Qian, I. Z. Wilson, A. Andreev, C.-C. Chen, and J.-H. Chu, Evidence for a strain-tuned topological phase transition in ZrTe5, Science Advances 5, eaav9771 (2019).
  * Fukui and Hatsugai [2007] T. Fukui and Y. Hatsugai, Quantum spin Hall effect in three dimensional materials: Lattice computation of $\mathbb{Z}_{2}$ topological invariants and its application to Bi and Sb, Journal of the Physical Society of Japan 76, 053702 (2007).
  * Asbóth  _et al._ [2016] J. K. Asbóth, L. Oroszlány, and A. Pályi,  _A Short Course on Topological Insulators: Band Structure and Edge States in One and Two Dimensions_, 1st ed., Lecture Notes in Physics, Vol. 909 (Springer International Publishing, BerlinXXX, 2016).
  * Li  _et al._ [2016a] X.-B. Li, W.-K. Huang, Y.-Y. Lv, K.-W. Zhang, C.-L. Yang, B.-B. Zhang, Y. Chen, S.-H. Yao, J. Zhou, M.-H. Lu, _et al._ , Experimental observation of topological edge states at the surface step edge of the topological insulator ZrTe5, Physical Review Letters 116, 176803 (2016a).
  * Manzoni  _et al._ [2016] G. Manzoni, L. Gragnaniello, G. Autès, T. Kuhn, A. Sterzi, F. Cilento, M. Zacchigna, V. Enenkel, I. Vobornik, L. Barba, _et al._ , Evidence for a strong topological insulator phase in ZrTe5, Physical Review Letters 117, 237601 (2016).
  * Wu  _et al._ [2016] R. Wu, J.-Z. Ma, S.-M. Nie, L.-X. Zhao, X. Huang, J.-X. Yin, B.-B. Fu, P. Richard, G.-F. Chen, Z. Fang, _et al._ , Evidence for topological edge states in a large energy gap near the step edges on the surface of zrte 5, Physical Review X 6, 021017 (2016).
  * Xiong  _et al._ [2017] H. Xiong, J. Sobota, S.-L. Yang, H. Soifer, A. Gauthier, M.-H. Lu, Y.-Y. Lv, S.-H. Yao, D. Lu, M. Hashimoto, _et al._ , Three-dimensional nature of the band structure of ZrTe5 measured by high-momentum-resolution photoemission spectroscopy, Physical Review B 95, 195119 (2017).
  * Lv  _et al._ [2018] Y.-Y. Lv, B.-B. Zhang, X. Li, K.-W. Zhang, X.-B. Li, S.-H. Yao, Y. Chen, J. Zhou, S.-T. Zhang, M.-H. Lu, _et al._ , Shubnikov–de haas oscillations in bulk zrt e 5 single crystals: Evidence for a weak topological insulator, Physical Review B 97, 115137 (2018).
  * Facio  _et al._ [2022] J. I. Facio, E. Nocerino, I. C. Fulga, R. Wawrzynczak, J. Brown, G. Gu, Q. Li, M. Mansson, Y. Sassa, O. Ivashko, _et al._ , Engineering a pure dirac regime in ZrTe5, arXiv preprint arXiv:2206.13957 10.48550/arXiv.2206.13957 (2022).
  * Tajkov  _et al._ [2022] Z. Tajkov, D. Nagy, K. Kandrai, J. Koltai, L. Oroszlány, P. Süle, Z. E. Horváth, P. Vancsó, L. Tapasztó, and P. Nemes-Incze, Revealing the topological phase diagram of ZrTe5 using the complex strain fields of microbubbles, npj Computational Materials 8, 177 (2022).
  * Stillwell  _et al._ [1989] E. Stillwell, A. Ehrlich, G. Kamm, and D. Gillespie, Effect of elastic tension on the electrical resistance of HfTe5 and ZrTe5, Physical Review B 39, 1626 (1989).
  * Zhang  _et al._ [2017] Y. Zhang, C. Wang, L. Yu, G. Liu, A. Liang, J. Huang, S. Nie, X. Sun, Y. Zhang, B. Shen, _et al._ , Electronic evidence of temperature-induced lifshitz transition and topological nature in zrte5, Nature communications 8, 1 (2017).
  * Niu  _et al._ [2017] J. Niu, J. Wang, Z. He, C. Zhang, X. Li, T. Cai, X. Ma, S. Jia, D. Yu, and X. Wu, Electrical transport in nanothick ZrTe5 sheets: From three to two dimensions, Physical Review B 95, 035420 (2017).
  * Li  _et al._ [2018] P. Li, C. Zhang, J. Zhang, Y. Wen, and X. Zhang, Giant planar hall effect in the Dirac semimetal ZrTe5-$\delta$, Physical Review B 98, 121108 (2018).
  * Liang  _et al._ [2018] T. Liang, J. Lin, Q. Gibson, S. Kushwaha, M. Liu, W. Wang, H. Xiong, J. A. Sobota, M. Hashimoto, P. S. Kirchmann, _et al._ , Anomalous hall effect in ZrTe5, Nature Physics 14, 451 (2018).
  * Li  _et al._ [2016b] Q. Li, D. E. Kharzeev, C. Zhang, Y. Huang, I. Pletikosić, A. Fedorov, R. Zhong, J. Schneeloch, G. Gu, and T. Valla, Chiral magnetic effect in ZrTe5, Nature Physics 12, 550 (2016b).
  * Liu  _et al._ [2016] Y. Liu, X. Yuan, C. Zhang, Z. Jin, A. Narayan, C. Luo, Z. Chen, L. Yang, J. Zou, X. Wu, _et al._ , Zeeman splitting and dynamical mass generation in dirac semimetal ZrTe5, Nature Communications 7, 1 (2016).
  * Chen  _et al._ [2015] R. Chen, S. Zhang, J. Schneeloch, C. Zhang, Q. Li, G. Gu, and N. Wang, Optical spectroscopy study of the three-dimensional Dirac semimetal ZrTe5, Physical Review B 92, 075107 (2015).
  * Zhang  _et al._ [2021] P. Zhang, R. Noguchi, K. Kuroda, C. Lin, K. Kawaguchi, K. Yaji, A. Harasawa, M. Lippmaa, S. Nie, H. Weng, _et al._ , Observation and control of the weak topological insulator state in zrte5, Nature communications 12, 1 (2021).
  * Xiang  _et al._ [2016] H. Xiang, B. Xu, J. Liu, Y. Xia, H. Lu, J. Yin, and Z. Liu, Quantum spin Hall insulator phase in monolayer WTe2 by uniaxial strain, AIP Advances 6, 095005 (2016).
  * Monserrat and Narayan [2019] B. Monserrat and A. Narayan, Unraveling the topology of ZrTe5 by changing temperature, Physical Review Research 1, 033181 (2019).
  * Tajkov  _et al._ [2019] Z. Tajkov, D. Visontai, L. Oroszlány, and J. Koltai, Uniaxial strain induced topological phase transition in bismuth–tellurohalide–graphene heterostructures, Nanoscale 11, 12704 (2019).
  * Tajkov  _et al._ [2020] Z. Tajkov, J. Koltai, J. Cserti, and L. Oroszlány, Competition of topological and topologically trivial phases in patterned graphene based heterostructures, Physical Review B 101, 235146 (2020).
  * Zhang  _et al._ [2022] C. Zhang, H. Pi, L. Zhou, S. Li, J. Zhou, A. Du, and H. Weng, Switchable topological phase transition and nonlinear optical properties in a ReC2H monolayer, Physical Review B 105, 245108 (2022).
  * Tang  _et al._ [2021] F. Tang, P. Wang, M. He, M. Isobe, G. Gu, Q. Li, L. Zhang, and J. H. Smet, Two-dimensional quantum Hall effect and zero energy state in few-layer ZrTe5, Nano letters 21, 5998 (2021).
  * Belke  _et al._ [2021] C. Belke, S. Locmelis, L. Thole, H. Schmidt, P. Behrens, and R. J. Haug, Thickness-dependent gap energies in thin layers of hfte5, 2D Mater. 8, 035029 (2021).
  * Zhuo  _et al._ [2022] W. Z. Zhuo, B. Lei, C. S. Zhu, Z. L. Sun, J. H. Cui, W. X. Wang, Z. Y. Wang, T. Wu, J. J. Ying, Z. J. Xiang, and X. H. Chen, Thickness-dependent electronic structure in layered ZrTe5 down to the two-dimensional limit, Phys. Rev. B 106, 085428 (2022).
  * Wang  _et al._ [2020] J. Wang, J. Niu, X. Li, X. Ma, Y. Yao, and X. Wu, Facile and fast growth of high mobility nanoribbons of ZrTe5, Chinese Physics B 29, 068102 (2020).
  * Weng  _et al._ [2014] H. Weng, X. Dai, and Z. Fang, Transition-metal Pentatelluride ZrTe5 and HfTe5: A paradigm for large-gap quantum spin Hall insulators, Physical Review X 4, 011002 (2014).
  * Magda  _et al._ [2015] G. Z. Magda, J. Pető, G. Dobrik, C. Hwang, L. P. Biró, and L. Tapasztó, Exfoliation of large-area transition metal chalcogenide single layers, Scientific Reports 5, 14714 (2015).
  * Velicky M  _et al._ [2018] G. E. Velicky M, Donnelly, W. R. Hendren, S. McFarland, D. Scullion, W. J. I. DeBenedetti, G. C. Correa, Y. Han, A. J. Wain, M. A. Hines, D. A. Muller, K. S. Novoselov, H. D. Abruna, R. M. Bowman, E. J. G. Santos, and F. Huang, Mechanism of Gold-Assisted exfoliation of Centimeter-Sized Transition-Metal dichalcogenide monolayers, ACS Nano 12, 10463 (2018).
  * Artacho  _et al._ [2008] E. Artacho, E. Anglada, O. Diéguez, J. D. Gale, A. García, J. Junquera, R. M. Martin, P. Ordejón, J. M. Pruneda, D. Sánchez-Portal, _et al._ , The SIESTA method; developments and applicability, Journal of Physics: Condensed Matter 20, 064208 (2008).
  * Soler  _et al._ [2002] J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, The SIESTA method for ab initio order-N materials simulation, Journal of Physics: Condensed Matter 14, 2745 (2002).
  * García  _et al._ [2020] A. García, N. Papior, A. Akhtar, E. Artacho, V. Blum, E. Bosoni, P. Brandimarte, M. Brandbyge, J. I. Cerdá, F. Corsetti, _et al._ , Siesta: Recent developments and applications, The Journal of chemical physics 152, 204108 (2020).
  * Fernández-Seivane  _et al._ [2006] L. Fernández-Seivane, M. A. Oliveira, S. Sanvito, and J. Ferrer, On-site approximation for spin–orbit coupling in linear combination of atomic orbitals density functional methods, Journal of Physics: Condensed Matter 18, 7999 (2006).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Physical Review B 47, 558 (1993).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical Review B 54, 11169 (1996).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Physical Review Letters 77, 3865 (1996).
  * Rivero  _et al._ [2015] P. Rivero, V. M. García-Suárez, D. Pereñiguez, K. Utt, Y. Yang, L. Bellaiche, K. Park, J. Ferrer, and S. Barraza-Lopez, Systematic pseudopotentials from reference eigenvalue sets for DFT calculations, Computational Materials Science 98, 372 (2015).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Physical Review B 13, 5188 (1976).
  * Papior [2021] N. Papior, sisl: v0.11.0 (2021).
  * Nagy [2021] D. Nagy, get_z2 (2021).
  * Klimeš  _et al._ [2011] J. Klimeš, D. R. Bowler, and A. Michaelides, Van der Waals density functionals applied to solids, Physical Review B 83, 195131 (2011).


# A Robust Alternative for Graph Convolutional Neural Networks via Graph
Neighborhood Filters

Victor M. Tenorio∗, Samuel Rey∗, Fer nando Gama†, Santiago Segarra†, and
Antonio G. Marques∗ Work in this paper was supported by the Spanish Grants
SPGRAPH (PID2019-105032GB-I00), FPU17/04520, and CAM PEJ-2020-AI/TIC-18964,
and the USA NSF award CCF-2008555. ∗Dept. of Signal Theory and Comms., King
Juan Carlos University, Madrid, Spain. † Dept. of ECE, Rice University,
Houston, USA. Contact author: antonio.garcia.marques(AT)urjc.es.

###### Abstract

Graph convolutional neural networks (GCNNs) are popular deep learning
architectures that, upon replacing regular convolutions with graph filters
(GFs), generalize CNNs to irregular domains. However, classical GFs are prone
to numerical errors since they consist of high-order polynomials. This problem
is aggravated when several filters are applied in cascade, limiting the
practical depth of GCNNs. To tackle this issue, we present the neighborhood
graph filters (NGFs), a family of GFs that replaces the powers of the graph
shift operator with $k$-hop neighborhood adjacency matrices. NGFs help to
alleviate the numerical issues of traditional GFs, allow for the design of
deeper GCNNs, and enhance the robustness to errors in the topology of the
graph. To illustrate the advantage over traditional GFs in practical
applications, we use NGFs in the design of deep neighborhood GCNNs to solve
graph signal denoising and node classification problems over both synthetic
and real-world data.

###### Index Terms:

Graph Neural Networks, Graph Filters, Robust GSP, Non-Euclidean Data,
Geometric Deep Learning.

##  I Introduction

The increasing complexity of current datasets, which oftentimes exhibit an
underlying irregular structure, compels us to develop new models capable of
learning efficiently from the observed data. One alternative to exploit the
irregular structure of the data at hand is provided by graph signal processing
(GSP) [1, 2, 3, 4, 5]. GSP is a rapidly growing field that assumes that the
properties of the data are closely related to its underlying irregular
structure, which can be accurately represented by a graph. A prominent family
of architectures provided by this discipline are graph convolutional neural
networks (GCNNs) [6, 7, 8], which leverage the generalization of the
convolution to irregular domains. In recent years, these architectures have
achieved state-of-the-art (SoA) performance in a wide range of applications
involving graph-structured data [9, 10, 11, 12, 13, 14, 15], including graph-
signal denoising, which is used as an illustrative test case in the manuscript
[16, 17, 18].

One of the key elements behind the success of GCNNs are graph filters (GFs)
[19, 20, 21], which are linear operators that employ the structure of the
graph to generalize the notion of classical convolution to graph signals. To
that end, GFs are defined as polynomials of the graph-shift operator (GSO), a
matrix encoding the topology of the observed graph. However, despite the
success of GCNNs and the well-known benefits associated with classical GFs,
the polynomial definition of the filters also comes with some limitations.
First, classical GFs are prone to numerical errors, a problem that is
aggravated when several filters are placed in cascade [22]. This numerical
instability can limit the design of deeper architectures. Another relevant
problem arises when there is uncertainty about the topology of the graph.
Because of their polynomial nature, GFs are sensitive to imperfections in the
observed graph, harming the performance of the subsequent GSP tasks [23, 24,
25, 26].

Motivated by the previous discussion, we present a new type of linear graph-
signal operators, referred to as neighborhood graph filters (NGFs), that
replace the powers of the GSO with $k$-hop adjacency matrices. These matrices
encode the topological information of $k$-hop neighborhoods by capturing the
existence of one or more shortest paths of a specific length between the nodes
of a given graph. As a result, the output signal generated by an NGF can be
interpreted as a linear combination of multiple signals, each of them
consisting in the aggregation of the input values at the nodes located at a
particular distance. We discuss the main properties of the proposed NGFs with
an especial focus on the numerical stability and the robustness to topology
perturbations. Furthermore, we exploit these filters to provide an alternative
design of GCNN that is employed in the unsupervised task of graph signal
denoising, and in the supervised problem of graph classification. The
performance of the resulting architecture is evaluated using both synthetic
and real-world datasets.

The remainder of the paper is organized as follows. Section II provides basic
concepts about GSP. Section III introduces the definition of NGFs and
discusses their properties. Section IV details the architecture resulting from
combining GCNNs with NGFs, and Section V provides numerical validation of the
presented architectures. Concluding remarks in Section VI wrap up the paper.

##  II Fundamentals of graph signal processing

This section introduces notation and reviews basic GSP concepts that are
leveraged throughout the paper.

Graphs. Let ${\mathcal{G}}=\\{{\mathcal{V}},{\mathcal{E}}\\}$ be a graph with
$N$ nodes collected in the set ${\mathcal{V}}$, and a set of edges
${\mathcal{E}}$ such that $(i,j)$ belongs to ${\mathcal{E}}$ if the nodes $i$
and $j$ are connected. For any given ${\mathcal{G}}$, its adjacency matrix is
represented by the ${N\times N}$ matrix ${\mathbf{A}}$ with nonzero elements
$[{\mathbf{A}}]_{ij}$ if and only if $(i,j)\in{\mathcal{E}}$. In other words,
the adjacency matrix encodes the 1-hop neighborhoods of the graph. Note that
for unweighted graphs, the entries of ${\mathbf{A}}$ are either $0$ or $1$.
Finally, we use ${\mathcal{N}}_{i}:=\\{j:(i,j)\in{\mathcal{E}}\\}$ to denote
the neighborhood of the node $i$ (i.e., the set of nodes that are linked to
$i$); $d_{i,j}$ to denote the minimum number of hops between nodes $i$ and
$j$; and $D=\max_{(i,j)}d_{i,j}$ to denote the diameter of ${\mathcal{G}}$,
which represents the length of the maximum shortest path present in the graph.

Graph signals. Signals observed on top of the graph are known as graph
signals. More formally, a graph signal is a function
$f:{\mathcal{V}}\rightarrow{\mathbb{R}}$ that can be represented as the vector
$[x_{1},...,x_{N}]^{T}\in{\mathbb{R}}^{N}$, where $x_{i}=f(i)$ denotes the
value of the signal observed at node $i$. When modeling and processing graph
signals, a key concept is the GSO ${\mathbf{S}}\in{\mathbb{R}}^{N\times N}$
[19], which is a matrix whose entries $(i,j)$ can be nonzero if and only if
$(i,j)\in{\mathcal{E}}$ or if $i=j$. The matrix ${\mathbf{S}}$ captures the
topology of the graph but makes no assumptions about the nonzero entries.
Typical choices for the GSO are the adjacency matrix ${\mathbf{A}}$ [19], the
graph Laplacian
${\mathbf{L}}=\text{diag}({\mathbf{A}}\mathbf{1})-{\mathbf{A}}$ [1], and their
respective generalizations. The GSO ${\mathbf{S}}$ represents a linear
transformation that can be computed locally at each node by aggregating the
values of the input signal within the one-hop neighborhood of each of the
nodes.

GFs. One of the most important tools in the context of GSP are GFs [19, 20].
GFs are linear graph signal operators that are defined as polynomials of the
GSO. Mathematically, given an $N\times N$ matrix ${\mathbf{H}}$, the operator
represented by the application of ${\mathbf{H}}$ to a graph signal is a GF if
the matrix can be written as

${\mathbf{H}}=\sum_{k=0}^{K-1}h_{k}{\mathbf{S}}^{k},$ (1)

where $K-1$ denotes the degree of the filter and
${\mathbf{h}}=[h_{0},...,h_{K-1}]$ collects the values of each of the $K$
filter coefficients. Since each application of ${\mathbf{S}}$ entails only
exchanges among one-hop neighbors, when applying ${\mathbf{S}}^{k}$ to an
input ${\mathbf{x}}$, the signal is being diffused across one-hop neighbors
$k$ times. Effectively, this leads to a diffusion over a $k$-hop neighborhood
[20]. This readily implies that the output of a generic GF ${\mathbf{H}}$ to
an input ${\mathbf{x}}$ can be written as
${\mathbf{y}}={\mathbf{H}}{\mathbf{x}}=\sum_{k=0}^{K}h_{k}({\mathbf{S}}^{k}{\mathbf{x}})=h_{0}{\mathbf{x}}+h_{1}{\mathbf{S}}{\mathbf{x}}+h_{2}{\mathbf{S}}^{2}{\mathbf{x}}+...$,
i.e., a linear combination of graph signals, each of them corresponding to the
original input diffused across neighborhoods of increasing size [20]. The
definition in (1) is easy to understand in the vertex domain, can be related
to the classical definition of convolution and linear time-invariant (LTI)
systems, has a neat spectral interpretation (see, e.g., [19] for the
definition of the graph Fourier transform for signals and filters), and has
been effectively used in a number of problems (such as denoising,
deconvolution, or signal reconstruction, to name a few). However, one of its
main limitations arises when $K$ is large, since high powers of ${\mathbf{S}}$
may render the filter (or the output) numerically ill-defined [20, 22]. This
issue is not present in classical LTI systems because the associated shift in
the time domain is isometric, so that the energy (norm) of a signal is
preserved regardless of the number of times the shift is applied. In the
following sections, we explore some of these issues in further detail.

##  III Linear neighborhood graph filters

This section presents NGFs, a new type of linear operator for graph signals.
The motivation for NGFs is to preserve most of the intuition present in
classical polynomial GFs, while bypassing some of the numerical problems
associated with high-degree filters.

Let us set ${\mathbf{S}}={\mathbf{A}}$ and suppose that ${\mathcal{G}}$ is
unweighted, so that the entries of ${\mathbf{A}}$ are binary. Then, it is
well-known that ${\mathbf{A}}^{k}$ encodes the number of $k$-hop paths between
any pair of nodes. For example, if $[{\mathbf{A}}^{2}]_{i,j}=3$, then, there
are 3 paths of length 2 connecting nodes $i$ and $j$. Hence, for unweighted
graphs, it follows that the application of ${\mathbf{A}}^{k}$ mixes the
information within nodes that are at most $k$ hops away. It also demonstrates
that as $k$ increases, the entries of ${\mathbf{A}}^{k}$ grow very large (the
number of paths increases exponentially) leading to numerical issues. Our
approach in this section is to replace ${\mathbf{A}}^{k}$ with a matrix that,
while preserving the notion of $k$-neighborhood, does not grow arbitrarily
large with $k$.

To be precise, let us start by defining the $k$-hop adjacency matrices
${\mathbf{A}}_{{\mathcal{N}}(k)}\in\\{0,1\\}$ as an $N\times N$ matrix whose
entry $[{\mathbf{A}}_{{\mathcal{N}}(k)}]_{ij}=1$ only if the nodes $i$ and $j$
are connected by at least one shortest path of exactly length $k$ (i.e., if
the distance between $i$ and $j$ satisfies $d_{i,j}=k$). Note that this
definition implies ${\mathbf{A}}_{{\mathcal{N}}(0)}={\mathbf{I}}$ and
${\mathbf{A}}_{{\mathcal{N}}(1)}={\mathbf{A}}$. With this notation at hand, we
can then define the NGF ${\mathbf{H}}_{\mathcal{N}}$ as a linear operator for
graph signals that can be written as [cf. (1)]

${\mathbf{H}}_{{\mathcal{N}}}=\sum_{k=0}^{K-1}h_{k}{\mathbf{A}}_{{\mathcal{N}}(k)},$
(2)

where ${\mathbf{h}}=[h_{0},...,h_{K-1}]$ are the filter coefficients. As done
in the previous section, when applied to an input graph signal ${\mathbf{x}}$,
the NGF generates the output
${\mathbf{y}}={\mathbf{H}}_{{\mathcal{N}}}{\mathbf{x}}$ where

${\mathbf{y}}=\sum_{k=0}^{K-1}h_{k}{\mathbf{A}}_{{\mathcal{N}}(k)}{\mathbf{x}}=h_{0}{\mathbf{x}}+h_{1}{\mathbf{A}}_{{\mathcal{N}}(1)}{\mathbf{x}}+h_{1}{\mathbf{A}}_{{\mathcal{N}}(2)}{\mathbf{x}}+...$
(3)

Unlike classical GFs, NGFs are not prone to numerical errors since the entries
of ${\mathbf{A}}_{{\mathcal{N}}(k)}$ are never larger than 1 and the norm
$\|{\mathbf{A}}_{{\mathcal{N}}(k)}\|$ need not increase with $k$. As a result,
NGFs are less sensitive to numerical instability. A particular interesting
property of NGFs is that ${\mathbf{A}}_{{\mathcal{N}}(k)}=\mathbf{0}_{N\times
N}$ for all $k>D$, where $D$ represents the diameter of ${\mathcal{G}}$,
limiting the number of active filter coefficients to $D+1$. While this fact
once again demonstrates the stability of NGFs, it also shows a potential loss
on the expressiveness of the proposed filters, especially when $D$ is small.
This issue can be addressed by incorporating more expressive filters, e.g.,
node variant GFs [20] or edge variant GFs [22]. Nevertheless, the robustness
to numerical issues derived from using the filters
${\mathbf{H}}_{{\mathcal{N}}}$ as proposed in (2) is expected to be more
significant when the graph presents a high diameter $D$, as it is illustrated
in Section V.

Another relevant property of ${\mathbf{H}}_{{\mathcal{N}}}$ is that, in
computing each entry of the output, the value of the input signal at each node
is considered at most once. Indeed, if
$[{\mathbf{A}}_{{\mathcal{N}}(k^{\prime})}]_{ij}=1$ then
$[{\mathbf{A}}_{{\mathcal{N}}(k)}]_{ij}=0$ for all $k\neq k^{\prime}$. Hence,
assuming that ${\mathcal{G}}$ is a connected graph and $K-1\geq D$, we have
that

$\sum_{k=0}^{K-1}{\mathbf{A}}_{{\mathcal{N}}(k)}=\mathbf{1}_{N\times N},$ (4)

where $\mathbf{1}_{N\times N}$ denotes the $N\times N$ matrix of all ones.
Moreover, let ${\bar{\mathbf{A}}}={\mathbf{A}}+{\mathbf{E}}$ represent a
perturbed adjacency matrix, with ${\mathbf{E}}$ denoting an error matrix that
randomly removes or adds edges from ${\mathbf{A}}$, and let
${\bar{\mathbf{A}}}_{{\mathcal{N}}(k)}$ be the $k$-hop adjacency matrices of
the perturbed graph associated with ${\bar{\mathbf{A}}}$. Since (4) holds for
both ${\mathbf{A}}_{{\mathcal{N}}(k)}$ and
${\bar{\mathbf{A}}}_{{\mathcal{N}}(k)}$, if the true and the perturbed graphs
are connected graphs, when one link is removed several shortest paths are
destroyed and the same number of shortest paths are created, and thus, for
$K-1\geq D$ we have that

$\sum_{k=0}^{K-1}\left({\bar{\mathbf{A}}}_{{\mathcal{N}}(k)}-{\mathbf{A}}_{{\mathcal{N}}(k)}\right)=\sum_{k=0}^{K-1}{\mathbf{E}}_{{\mathcal{N}}(k)}=\mathbf{0}_{N\times
N},$ (5)

where ${\mathbf{E}}_{{\mathcal{N}}(k)}$ represents the error induced by the
matrix ${\mathbf{E}}$ in the $k$-hop adjacency matrix
${\mathbf{A}}_{{\mathcal{N}}(k)}$. From (5), it can be observed that if all
the filter coefficients are constant (i.e., $h_{k}=h$ for all $k$), the error
between the perturbed and the true NGF is given by

$\|{\bar{\mathbf{H}}}_{{\mathcal{N}}}-{\mathbf{H}}_{{\mathcal{N}}}\|^{2}=\left\|h\sum_{k=0}^{K-1}\left({\bar{\mathbf{A}}}_{{\mathcal{N}}(k)}-{\mathbf{A}}_{{\mathcal{N}}(k)}\right)\right\|^{2}=0,$
(6)

where it can be seen that for this particular setting the NGF is impervious to
topology perturbations. GFs with constant coefficients appear, for example, in
some message-passing applications [27, 28].

The preliminary result presented in (6) provides mathematical support to the
idea that NGFs are more robust to topology imperfections than classical GFs.
The perturbed GF with ${\bar{\mathbf{A}}}$ as the GSO is given by
${\bar{\mathbf{H}}}=\sum_{k=0}^{K-1}h_{k}{\bar{\mathbf{A}}}^{k}=\sum_{k=0}^{K-1}h_{k}({\mathbf{A}}+{\mathbf{E}})^{k}$,
so intuitively, it can be observed that the powers of ${\bar{\mathbf{A}}}$
will amplify the errors on ${\mathbf{A}}$, increasing the discrepancies
between ${\mathbf{H}}$ and ${\bar{\mathbf{H}}}$ as $k$ increases [26]. On the
other hand, regarding NGFs, when a path of length $k$ is perturbed due to
imperfections in the observed adjacency matrix ${\bar{\mathbf{A}}}$, the
redundancy of paths existing in most graphs suggests that, as $k$ increases,
it is more likely that an alternative path of the same length will also be
available. Therefore, it is expected that the difference between the matrices
${\mathbf{A}}_{{\mathcal{N}}(k)}$ and ${\bar{\mathbf{A}}}_{{\mathcal{N}}(k)}$
decreases for higher values of $k$. This intuition is evaluated numerically in
Figure 1, where the normalized errors
${\|{\bar{\mathbf{H}}}-{\mathbf{H}}\|_{F}^{2}}/{\|{\mathbf{H}}\|_{F}^{2}}$ and
${\|{\bar{\mathbf{H}}}_{{\mathcal{N}}}-{\mathbf{H}}_{{\mathcal{N}}}\|_{F}^{2}}/{\|{\mathbf{H}}_{{\mathcal{N}}}\|_{F}^{2}}$
are depicted for classical GFs and NGFs, respectively. The comparison is
carried out in two different random graphs: Erdős Rényi and small-world graphs
[29]. In both cases, it can be seen that the error of classical GFs increases
rapidly as $K$ grows.

This observation notwithstanding, the facts that matrices
${\mathbf{A}}_{{\mathcal{N}}(k)}$ are not simultaneously diagonalizable and
that they represent the presence of shortest paths (which is a non-
differentiable operation) render the theoretical characterization of the
robustness of ${\mathbf{H}}_{{\mathcal{N}}}$ a challenging task that will be
addressed in future works.

Figure 1: Evolution of the normalized error for classical GFs and NGFs as the
number of filter taps $K$ increases for Erdős-Rényi and small-world random
graphs. The result is the mean error of 100 realizations.

##  IV NGF-based GCNN

In this section, we introduce a natural extension of GCNNs, where the
classical GFs are replaced by the novel NGFs.

A graph neural network is a parametric non-linear function

$f_{\boldsymbol{\Theta}}({\mathbf{Z}}|{\mathcal{G}}):~{}{\mathbb{R}}^{N\times
F^{(0)}}\rightarrow{\mathbb{R}}^{N\times F^{(L)}},$ (7)

that depends on the graph ${\mathcal{G}}$. The matrix ${\mathbf{Z}}$
represents the input of the architecture, the learnable weights are collected
in $\boldsymbol{\Theta}=\\{\boldsymbol{\Theta}^{(\ell)}\\}_{\ell=1}^{L}$, $L$
represents the number of layers, and $F^{(\ell)}$ denotes the number of
features at layer $\ell$. When designing a graph neural network there are
several alternatives to account for the topology of the graph. One common
approach used in specific implementations of GCNNs exploits the message
passing operation. The resulting architecture is given by the following
recursion

${\mathbf{X}}^{(\ell)}=\sigma\big{(}{\mathbf{A}}{\mathbf{X}}^{(\ell-1)}\boldsymbol{\Theta}^{(\ell)}\big{)},$
(8)

where $\sigma(\cdot)$ is an element-wise non-linear function typically known
as activation function, the weights $\boldsymbol{\Theta}^{(\ell)}$ are
matrices of dimension $F^{(\ell-1)}\times F^{(\ell)}$, and
${\mathbf{X}}^{(\ell)}\in{\mathbb{R}}^{N\times F^{(\ell)}}$ denotes the output
of the $\ell$-th layer with ${\mathbf{X}}^{(0)}={\mathbf{Z}}$ and
${\mathbf{X}}^{(L)}$ being the input and the output of the architecture,
respectively. Note that we can interpret the computation of
${\mathbf{X}}^{(\ell)}$ as first left-multiplying the input of the previous
layer by the adjacency matrix of the graph ${\mathbf{A}}$, which combines the
signal values of neighboring nodes, and then right-multiplying by the
learnable weights $\boldsymbol{\Theta}^{(\ell)}$, so the architecture learns
to combine different features from previous layers.

The main disadvantage of the architecture presented in (8) is that, at each
layer, the convolution is carried out only in the 1-hop neighborhoods. This
effect is mitigated by stacking $L$ layers, since it is roughly equivalent to
perform a convolution on the $L$-hop neighborhood. However, this alternative
produces an undesirable coupling between the depth of the GCNN and the size of
the neighborhood where the convolution is applied.

Figure 2: Comparison of the performance of the proposed NGF against a
classical GF setting in the synthetic graph signal denoising task. (a) Median
error as a function of training epochs for a bandlimited signal corrupted by
noise with power $0.1$. (b) Median error achieved by the architectures as a
function of noise power when denoising a diffused white signal (for each
experiment, the number of epochs corresponds to the one that minimizes the
error). In both cases, the signal is normalized to have unit norm.

One way to avoid the aforementioned coupling is to replace the matrix
${\mathbf{A}}$ with a graph filter ${\mathbf{H}}$, resulting in an
architecture implemented by the recursion

${\mathbf{X}}^{(\ell)}=\sigma\big{(}{\mathbf{H}}^{(\ell)}{\mathbf{X}}^{(\ell-1)}\boldsymbol{\Theta}^{(\ell)}\big{)}.$
(9)

Note that ${\mathbf{H}}^{(\ell)}$ is a classical filter of the form (1) and
the filter coefficients can be: 1) fixed in advance as a prior or 2) learned
by the architecture, giving rise to two different architectures. Although GFs
decouple the depth of the architecture from the range of the convolution, the
drawbacks introduced in Section III still impose some limitations on the
architecture. One evident restriction affects the depth of the architecture.
Stacking layers is similar to applying several GFs in cascade, so the
numerical issues derived from high-order polynomials can give rise to
exploding or vanishing gradients (especially when ReLUs are used as activation
functions). To circumvent this problem, in this work we replace the classical
GF with the NGF, so the neighborhood GCNN (NGCNN) is given by

${\mathbf{X}}^{(\ell)}=\sigma\big{(}{\mathbf{H}}_{{\mathcal{N}}}^{(\ell)}{\mathbf{X}}^{(\ell-1)}\boldsymbol{\Theta}^{(\ell)}\big{)},$
(10)

where ${\mathbf{H}}_{{\mathcal{N}}}^{(\ell)}$ is a filter of the form (2) and
the filter coefficients can be either fixed in advance or learned by the
architecture. The proposed NGCNN preserves the structure of the graph-aware
per-layer linear transformation, so the weights are only learning to mix the
different features of the input matrices while the relation between the signal
values at the different nodes is determined by the NGF. Note that, from a GSP
perspective, since the matrix ${\mathbf{X}}^{(\ell-1)}$ is interpreted as
$F^{(\ell-1)}$ different graph signals, the product
${\mathbf{H}}_{{\mathcal{N}}}^{(\ell)}{\mathbf{X}}^{(\ell-1)}$ is seen as a
node-domain convolution. Moreover, thanks to the filters
${\mathbf{H}}_{{\mathcal{N}}}^{(\ell)}$, the depth of the architecture and the
range of the convolutions are completely decoupled since the NGF are less
prone to numerical issues, and thus, they do not impose a limiting factor in
the depth of the architecture. On the other hand, the maximum degree of the
NGF at each layer is $D$, the diameter of the graph.

Next, we evaluate the performance of the proposed NGCNN in different settings
and numerically illustrate how NGFs mitigate some of the typical limitations
of GCNN.

##  V Numerical Simulations

We analyze the performance of the proposed NGF in two different graph-related
problems: signal denoising using synthetic data (Section V-A) and node
classification in three real-world citation networks (Section V-B). The python
code used to run the experiments, which describes in detail all the settings
of the proposed architecture, is available
online111https://github.com/vmtenorio/NeighborhoodGF.

###  V-A Graph signal denoising

The goal in graph-signal denoising is to recover the original graph signal
${\mathbf{x}}\in{\mathbb{R}}^{N}$ given the noisy graph-signal observation
${\mathbf{y}}={\mathbf{x}}+{\mathbf{w}}$, with
${\mathbf{w}}\in{\mathbb{R}}^{N}$ representing a noise vector. To that end, we
approach the denoising problem as in [17, 30] by minimizing

$\displaystyle\\!\hat{\boldsymbol{\Theta}}=\operatornamewithlimits{argmin}_{\boldsymbol{\Theta}}$
$\displaystyle\|{\mathbf{y}}-{\mathbf{f}}_{\boldsymbol{\Theta}}({\mathbf{Z}}|{\mathcal{G}})\|_{2}^{2},$
(11)

where the entries of the input matrix ${\mathbf{Z}}$ are randomly sampled from
a zero-mean unit-variance normal distribution. The problem (11) is minimized
by running gradient descent for a fixed number of epochs. After estimating the
weights $\hat{\boldsymbol{\Theta}}$, the denoised signal is given by
${\hat{\mathbf{x}}}={\mathbf{f}}_{\hat{\boldsymbol{\Theta}}}({\mathbf{Z}}|{\mathcal{G}})$.
The motivation behind this approach is that the proposed architecture is
capable of learning the signal ${\mathbf{x}}$ faster than the noise
${\mathbf{w}}$ and, thus, early stopping can be applied to learn most of the
signal without learning too much noise.

To analyze the performance of the proposed filters, we created synthetic graph
signals defined over a random stochastic-block-model (SBM) [29] graph. In this
graph-generative model, nodes are divided into communities and edges are
randomly drawn between nodes of the same community independently with a
probability of 0.3, and between nodes of different communities, also
independently, with a probability of 0.0075. The graph contains 256 nodes and
8 or 4 equally-sized communities for test cases 1 and 2, respectively.

Figure 3: Comparison of the performance of the proposed NGF against a
classical GF setting in three real-world scenarios. (a) Test classification
accuracy for the three considered datasets as a function of the number of
filter coefficients used in the network. (b) Performance of the architectures
using the Citeseer dataset, for different numbers of filter taps ($K$) as a
function of the amount of perturbations introduced to the graph.

Test case 1. The results for the first test case are shown in Figure 2. The
signal ${\mathbf{x}}$ is created by diffusing a zero-mean white signal
${\mathbf{b}}\in{\mathbb{R}}^{N}$ over the graph as
${\mathbf{x}}={\mathbf{H}}{\mathbf{b}}$ for the case “Input GF”, and
${\mathbf{x}}={\mathbf{H}}_{{\mathcal{N}}}{\mathbf{b}}$ for “Input NGF”. Both
filters ${\mathbf{H}}$ and ${\mathbf{H}}_{{\mathcal{N}}}$ have random
coefficients ${\mathbf{h}}$, drawn from a white uniform distribution between 0
and 1 and are normalized so that ${\mathbf{h}}^{T}\mathbf{1}=1$. We add zero-
mean white Gaussian noise ${\mathbf{w}}$ with normalized power of $0.1$, and
generate 200 different signals ${\mathbf{x}}$. The results shown in the figure
report the median error
$\|{\mathbf{x}}-{\hat{\mathbf{x}}}\|_{2}^{2}/\|{\mathbf{x}}\|_{2}^{2}$ across
the 200 realizations as the number of epochs increases for four different
2-layer architectures: one using the classical GF (“Arch. GF”), another one
using the NGF (“Arch. NGF”), and two SoA architectures: Graph Convolutional
Networks [10] (GCN) and Simple Graph Convolution [31] (SGC). It can be
observed that the error starts decreasing, reaches a minimum, and then
increases. This latter increase is due to the architecture starting to fit the
noise, so early stopping is needed to obtain the denoised
${\hat{\mathbf{x}}}$. It can also be seen that the performance is tightly
coupled to the signal generation method. If we generate ${\mathbf{x}}$ using
${\mathbf{H}}$ then the architecture with the best performance is the one
using a classical GF, beating both SoA architectures in terms of minimum
error. However, when generating the signal with the NGF, the architecture with
superior performance and the only one able to denoise the signal is the one
using an NGF. As a result, NGCNN is the only architecture to effectively
denoise the signals in the two scenarios considered.

Test case 2. The signal ${\mathbf{x}}$ is generated as in the previous test
case. The results are shown in Figure 2, where we depict the evolution of the
normalized error achieved by each architecture as the normalized power of the
noise increases. The figure reveals that NGFs and classic GFs behave similarly
when the input is generated using a classic GF, with GCN outperforming both
architectures for high noise power values. However, NGF clearly outperforms
classic GF in all cases when the signal is generated with a diffusion process
that utilizes an NGF, illustrating that NGFs are more flexible linear
operators. NGCNN also outperforms GCN in low noise power settings, while GCN
achieves a lower minimum error by a small margin when the noise power is
greater than 0.2

###  V-B Node Classification

We analyzed the performance of the proposed graph filters using 3 real-world
datasets, where the graph represents a network of publications (nodes are
published articles and edges represent citations between them) and the goal is
to classify each node into a category that denotes the topic of the
publication. The graph signals in each node indicate the presence or absence
of words from a dictionary. Further details of these datasets can be found in
Table I and [32, 33, 34].

To perform the classification task, the output of the architectures is given
by ${\mathbf{X}}^{(L)}\in{\mathbb{R}}^{N\times M}$, where $M$ denotes the
total number of classes. The non-linearity in the last layer is given by the
softmax function so that the features in each node are interpreted as a
measurement of the probability for the specific node to belong to a certain
class. The architectures are trained with the cross-entropy loss, a typical
loss function used for classification problems. To avoid numerical issues
related to the powers of the GSO, in these experiments we normalize the GSO by
its largest eigenvalue ${\tilde{\mathbf{S}}}={\mathbf{S}}/\lambda_{max}$.

TABLE I: Datasets used for node classification.

Dataset name | #Nodes | #Features | #Classes | Radius* | Diameter*
---|---|---|---|---|---
Cora [32] | 2708 | 1433 | 7 | 10 | 19
Citeseer [33] | 3327 | 3703 | 6 | 15 | 28
Pubmed [34] | 19717 | 500 | 3 | 10 | 18

* Metrics refer to the largest connected component of the graph.

Test case 3. With the configuration described above, the accuracy over the set
of test nodes obtained in each dataset for both types of GFs can be seen in
Figure 3. The image shows the performance of the GCNN and the NGCNN as the
order of the filters increases. It can be seen that, for $K=2$, the
performance of both filters is similar. This is not surprising since filters
of this order only take into account the 1-hop neighborhood, so in this case
${\mathbf{H}}={\mathbf{H}}_{{\mathcal{N}}}$. Nevertheless, as $K$ increases
the behavior aligns with the discussion provided in Section III. The
performance of the NGF improves while the performance of the classical GF
deteriorates or improves at a much slower rate. This illustrates how NGFs are
more robust to numerical issues related to higher-order filters This can be
seen especially clear for the Citeseer dataset, which has the graph with the
highest diameter among the three.

Test case 4. In this case the focus is on the Citeseer dataset, and we analyze
the performance of the network as we increase the perturbation introduced in
the graph; see Figure 3. The effect of the perturbation is measured for
different values of the number of filter coefficients $K$, as indicated in the
legend. The perturbations consist in randomly removing and creating links in
the original graph, as a percentage of the total number of existing links.
NGFs outperform classic GF in all the tested settings and the performance of
NGFs remains approximately constant independently of the perturbation
introduced, while the accuracy of the architectures with classical GFs
decreases. This result is aligned with the discussion presented throughout
this paper, showcasing that NGFs are more robust to errors in the topology of
the given graph and further motivating the development of associated
theoretical results.

##  VI Conclusions

This paper proposed graph neural network architectures for graph signals based
on NGFs, a new type of GFs whose definition leverages the (adjacency) matrices
encoding the $k$-hop neighborhood of the nodes of the graph. Compared with
classical GFs, the proposed NGFs architectures are i) more stable when the
number of filter coefficients increases and ii) more robust against graph
perturbations, rendering this new technique more suitable for problems where
the graph is not known with certainty. We validated these claims through
experimental results in synthetic datasets and, more conspicuously, in real-
world datasets with larger graphs in terms of both size and diameter. Future
work includes the theoretical characterization of the robustness of NGFs as
well as additional simulations in real-world tasks.

## References

  * [1] 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 Process. Mag., vol. 30, no. 3, pp. 83–98, May 2013\.
  * [2] E.D. Kolaczyk and G. Csárdi,  Statistical analysis of network data with R, vol. 65,  Springer, 2014.
  * [3] A.G. Marques, N. Kiyavash, J.M.F. Moura, D. Van De Ville, and R. Willett,  “Graph Signal Processing: Foundations and Emerging Directions [From the Guest Editors],”  IEEE Signal Process. Mag., vol. 37, no. 6, pp. 11–13, 2020.
  * [4] A. Ortega, P. Frossard, J. Kovačević, J.M.F. Moura, and P. Vandergheynst,  “Graph signal processing: Overview, challenges, and applications,”  Proc. IEEE, vol. 106, no. 5, pp. 808–828, 2018.
  * [5] P.M. Djurić and C. Richard,  Cooperative and Graph Signal Processing,  Academic Press, 2018.
  * [6] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini,  “The graph neural network model,”  IEEE Trans. Neural Netw., vol. 20, no. 1, pp. 61–80, Jan. 2009\.
  * [7] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst,  “Geometric deep learning: Going beyond euclidean data,”  IEEE Signal Process. Mag., vol. 34, no. 4, pp. 18–42, July 2017\.
  * [8] F. Gama, E. Isufi, G. Leus, and A. Ribeiro,  “Graphs, convolutions, and neural networks: From graph filters to graph neural networks,”  IEEE Signal Process. Mag., vol. 37, no. 6, pp. 128–138, Nov. 2020\.
  * [9] M. Defferrard, X. Bresson, and P. Vandergheynst,  “Convolutional neural networks on graphs with fast localized spectral filtering,”  in Conf. Neural Inform. Process. Syst., 2016, pp. 3844–3852.
  * [10] T.N. Kipf and M. Welling,  “Semi-Supervised Classification with Graph Convolutional Networks,”  in Int. Conf. Learn. Representations, Toulon, FR, 24–26 April, 2017\.
  * [11] W.L. Hamilton, Z. Ying, and J. Leskovec,  “Inductive representation learning on large graphs.,”  in Conf. Neural Inform. Process. Syst., 2017, pp. 1024–1034.
  * [12] F. Gama, A. G. Marques, G. Leus, and A. Ribeiro,  “Convolutional neural network architectures for signals supported on graphs,”  IEEE Trans. Signal Process., vol. 67, no. 4, pp. 1034–1049, Feb 2019.
  * [13] A. Chowdhury, G. Verma, C. Rao, A. Swami, and S. Segarra,  “Unfolding WMMSE using graph neural networks for efficient power allocation,”  IEEE Trans. Wireless Commun., vol. 20, no. 9, pp. 6004–6017, 2021\.
  * [14] T.M. Roddenberry, N. Glaze, and S. Segarra,  “Principled simplicial neural networks for trajectory prediction,”  in Intl. Conf. Machine Learn. (ICML), 2021.
  * [15] G. Cutura, B. Li, A. Swami, and S. Segarra,  “Deep demixing: Reconstructing the evolution of epidemics using graph neural networks,”  in European Signal Process. Conf. (EUSIPCO), Dublin, IR, 23–27 August, 2021.
  * [16] S. Chen, A. Sandryhaila, J.M.F. Moura, and J. Kovacevic,  “Signal denoising on graphs via graph filtering,”  in Global Conf. Signal and Info. Process. (GlobalSIP), 2014, pp. 872–876.
  * [17] S. Rey, A.G. Marques, and S. Segarra,  “An underparametrized deep decoder architecture for graph signals,”  in IEEE Intl. Wrksp. Computat. Advances Multi-Sensor Adaptive Process. (CAMSAP). IEEE, 2019, pp. 231–235.
  * [18] M. Onuki, S. Ono, M. Yamagishi, and Y. Tanaka,  “Graph signal denoising via trilateral filter on graph spectral domain,”  IEEE Trans. Signal Inf. Process. Netw., vol. 2, no. 2, pp. 137–148, 2016.
  * [19] A. Sandryhaila and J.M.F. Moura,  “Discrete signal processing on graphs,”  IEEE Trans. Signal Process., vol. 61, no. 7, pp. 1644–1656, Apr. 2013.
  * [20] S. Segarra, A.G. Marques, and A. Ribeiro,  “Optimal graph-filter design and applications to distributed linear network operators,”  IEEE Trans. Signal Process., vol. 65, no. 15, pp. 4117–4131, Aug 2017.
  * [21] F.M. Bianchi, D. Grattarola, L. Livi, and C. Alippi,  “Graph neural networks with convolutional arma filters,”  IEEE Trans. Pattern Anal. Mach. Intell., 2021.
  * [22] M. Coutino, E. Isufi, and G. Leus,  “Advances in distributed graph filtering,”  IEEE Trans. Signal Process., vol. 67, no. 9, pp. 2320–2333, 2019\.
  * [23] S. Segarra and A. Ribeiro,  “Stability and continuity of centrality measures in weighted graphs,”  IEEE Trans. Signal Process., vol. 64, no. 3, pp. 543–555, 2016\.
  * [24] E. Ceci and S. Barbarossa,  “Graph signal processing in the presence of topology uncertainties,”  IEEE Trans. Signal Process., vol. 68, pp. 1558–1573, 2020.
  * [25] A. Natali, M. Coutino, and G. Leus,  “Topology-aware joint graph filter and edge weight identification for network processes,”  in Intl Worksp. Machine Learn. Signal Process. (MLSP). IEEE, 2020, pp. 1–6.
  * [26] S. Rey and A.G. Marques,  “Robust graph-filter identification with graph denoising regularization,”  in IEEE Intl. Conf. on Acoustics, Speech and Signal Process., 2021, pp. 5300–5304.
  * [27] L. Zhang, D. Xu, A. Arnab, and P. H.S. Torr,  “Dynamic graph message passing networks,”  in Intl. Conf. Comput. Vision and Pattern Recognition, June 2020\.
  * [28] J. Dauwels,  “On variational message passing on factor graphs,”  in IEEE Intl. Symp. Info. Theory, 2007, pp. 2546–2550.
  * [29] M. Newman,  Networks,  Oxford University Press, 2018.
  * [30] T.H. Do, D.M. Nguyen, and N. Deligiannis,  “Graph auto-encoder for graph signal denoising,”  in IEEE Intl. Conf. on Acoustics, Speech and Signal Process. IEEE, 2020, pp. 3322–3326.
  * [31] F. Wu, A. Souza, T. Zhang, C. Fifty, T. Yu, and K. Weinberger,  “Simplifying graph convolutional networks,”  in Intl. Conf. Machine Learn. (ICML), 09–15 Jun 2019, vol. 97, pp. 6861–6871.
  * [32] A. Mccallum, K. Nigam, and J. Rennie,  “Automating the Construction of Internet Portals,”  Info. Retrieval, Springer, p. 127–163, 03 2000.
  * [33] C.L. Giles, K.D. Bollacker, and S. Lawrence,  “Citeseer: An automatic citation indexing system,”  in ACM Conf. Digital Libraries, 1998, p. 89–98.
  * [34] P. Sen, G. Namata, M. Bilgic, L. Getoor, B. Galligher, and T. Eliassi-Rad,  “Collective classification in network data,”  AI Mag., vol. 29, no. 3, pp. 93, Sep. 2008.


# HUT: Enabling High-UTility, Batched Queries
under Differential Privacy Protection for Internet-of-Vehicles

###### Abstract

The emerging trends of Internet-of-Vehicles (IoV) demand centralized servers
to collect/process sensitive data with limited computational resources on a
single vehicle. Such centralizations of sensitive data demand practical
privacy protections. One widely-applied paradigm, Differential Privacy, can
provide strong guarantees over sensitive data by adding noises. However,
directly applying DP for IoV incurs significant challenges for data utility
and effective protection. We observe that the key issue about DP-enabled
protection in IoV lies in how to synergistically combine DP with special
characteristics of IoV, whose query sequences are usually formed as unbalanced
batches due to frequent interactions between centralized servers and edge
vehicles.

To this end, we propose HUT, a new algorithm to enable High UTility for DP-
enabled protection in IoV. Our key insight is to leverage the inherent
characteristics in IoV: the unbalanced batches. Our key idea is to aggregate
local batches and apply Order Constraints, so that information loss from DP
protection can be mitigated. We evaluate the effectiveness of HUT against the
state-of-the-art DP protection mechanisms. The results show that HUT can
provide much lower information loss by 95.69% and simultaneously enable strong
mathematically-guaranteed protection over sensitive data.

##  1 Introduction

Internet-of-Vehicles (IoV) is considered a promising solution to resolve the
limited computational power of a single vehicle. IoV centralizes statistics
from edge vehicles in a power server for computation. However, the
centralizations of these statistics may contain sensitive data, and
substantially requires privacy protections in practice [1, 2, 3, 4, 5, 6].
Differential Privacy (DP) [7], as the dominant mechanism for privacy-
preserving computation, provides strong guarantees over sensitive data by
statistically adding noises. Though DP is a promising solution for the privacy
concerns in IoV, directly applying DP for IoV incurs significant challenges in
terms of both data utility and effective protection.

We address the key issue behind applying DP for IoV: how to synergistically
combine DP with special characteristics of IoV, whose query sequences are
usually formed as unbalanced batches in practice. The unique characteristic of
IoV is the unbalanced sequence of requests (denoted as batches). These batches
are formed due to a practical issue: vehicles require frequent queries of
relevant data for in-time decision-makings [8]. Hence, how to strike a balance
between privacy protection and data utility is a major challenge for applying
DP for IoV [9]. This is because privacy over-protection can cause low data
utility.

Prior works aim to balance the tradeoffs between data utility and privacy
protection, but fail to address the issues of batches within IoV. [10]
introduces Constraint Inference, which applies a constraint on the order of a
sorted data set to decrease unnecessary noise; [11] introduces Micro-
Aggregation, which uses a fixed-size cluster to cluster data set before adding
noise, to reduce the information loss; and [12] introduces K-aggregation,
which aggregates small-value data to form new data when the summation of
small-value data exceeds a certain threshold. The ignorance prohibits DP to be
applied effectively in the context of IoV, which demands new insights for
bridging DP with IoV.

We propose HUT, a new algorithm to enable High UTility for DP-enabled
protection in the context of IoV. Our key insight is to take advantage of the
fact that, the interactions between centralized servers and edge vehicles are
usually frequent and fine-grained, which are combined into a series of
unbalanced batches. To exploit the characteristic, our key idea is to
aggregate local batches and apply Order Constraints, so that information loss
from DP protection can be mitigated.

We quantitatively evaluate the effectiveness of HUT against the state-of-the-
art DP protection mechanism. Our results show that HUT can provide much lower
information loss by 95.69% and simultaneously enable strong mathematically-
guaranteed protection of sensitive data. Meanwhile, HUT does not influence the
strong guarantee of DP design, which ensures effective protection of sensitive
data from users.

##  2 Design Overview

We give an overview of our HUT design for ensuring DP protection with high
data utility. HUT consists of three steps: 1) we leverage Micro-aggregation
and K-means algorithm to group small-value batches together so that we can
mitigate the negative effects of uniform DP protection; 2) we deploy
Differential Privacy (DP) to preserve the privacy of the pre-processed data;
3) we exploit Order Constraint (OC) to re-sort DP-protected data so that we
can achieve high data utility. Figure 1 compares HUT with two other state-of-
the-art DP protection methods, and we elaborate the differences in details as
follows.

Fig. 1: An illustration of three DP-protection methods. ➊ refers to HUT (our
proposal), in which the details of the noisy-data processing method(OC)[10] is
illustrated; Two state-of-the-art methods are illustrated (➋ [11] and ➌ [12]).

1) Micro-Aggregation: HUT first applies Micro-Aggregation to cluster small-
value raw data, so that we can mitigate the negative impacts from DP in terms
of data utility. This is because DP adds uniform random noises to the dataset,
and small-value data will have a higher Signal-to-Noise Ratio (SNR), compared
with large-value data. Therefore, the data utility is impacted dependently
according to its exact values. Therefore, a mechanism to mitigate such impacts
is essential. We achieve so as follows. Before directly using DP techniques on
raw data for protection, HUT sets a threshold to divide small-value data and
the rest into different batches (as To-be-clustered batches and Not-clustered
batches). Next, we apply K-Mean clustering on data within the threshold (the
To-be-clustered batches), to automatically group such small-value data further
into smaller intra-batches (i.e. Micro-Aggregation).

Compared with [11], HUT effectively reduces the sensitivity and variations of
small-value data as they are upper bounded by the centroid of the respective
cluster. In [11] (shown as ➋ in Figure 1), a fixed-size clustering method is
applied and this results in a consistent distance function on a total order
relation but incurs large bias on the dataset which has highly uneven data
distribution. Such bias on unevenly distributed data could harm the integrity
of the dataset, leading to additional information loss. Thus, on the contrary,
we leverage the K-means clustering algorithm to perform the clustering to
mitigate the bias.

2) DP Protection: we then add Laplace noise on the data to achieve
$\epsilon$-Differential Privacy. DP aims to minimize the changes in query
answers when a record is changed in the data set by adding a randomized
function (subject to $\epsilon$-differential privacy on the data). Hereby, we
elaborate $\epsilon$-Differential Privacy via a mathematical description of a
randomized function $N_{f}$ (shown below).

###### Theorem 1

(Randomized Function subject to $\epsilon$-DP) We regard a randomized function
$N_{f}$ as satisfying $\epsilon$-DP, assuming that $\forall$ data set $X_{1}$
and $X_{2}$ with at most one element difference, and $\forall S\subseteq
Range(K_{f})$, the $N_{f}$ function obey the below in-equation:

$Pr[N_{f}(X_{1})\in S]\leq exp(\epsilon)\times Pr[N_{f}(X_{2})\in S]$

The above-randomized function is commonly realized by adding Laplace
Distributed Noise $\epsilon$, which is inversely proportional to the noise
magnitude (If $\epsilon$ is smaller, the magnitude of the noise is larger). In
HUT, we try out various $\epsilon$ settings to testify the proposed method
feasibility in different protection scales.

3) Order Constraint: HUT finally utilizes a noisy-data processing method,
Order Constraint(OC), to minimize unnecessary noise. A general feature of
noisy data is that after adding noise(especially large-scale noise) on a
previously sorted dataset, the noisy dataset often turns out to be disordered.
[10] introduces OC to re-process the noisy data to ensure the same order as
before. Specifically, the Min-Max formulas (constraint inference) [13] is used
to remove the unnecessary noise, applying least squares regression to change
the data values in the noisy dataset to make sure it still follows the
previous sorted order. The removed noise is proved to be unnecessary and does
not influence the effect of dataset privacy protection but can have a
significant improvement on reducing unnecessary information loss. Thus, we
apply this method under the context of DP protection as a step to further
reducing unnecessarily-added DP noise while having no negative impacts on
privacy protection itself.

Novelty of HUT: HUT has three novelties, compared with the state-of-the-art DP
protection mechanisms. First, HUT applies two-stage Micro-aggregation to
mitigate the negative impacts on fine-grained batches, which is considered as
a unique characteristic of IoV; Second, HUT leverages K-Means Clustering,
rather than fixed-size clustering, to improve the utility of protected data;
and third, HUT selectively enables lightweight constraints to be synergistic
with other designs, for both DP protection and data utility .

Are there any negative impacts on DP protection: In short, there are no
negative impacts on DP protection when applying HUT: since DP protection and
K-means clustering are not in conflict with the mathematical guarantee of DP,
we focus on the Micro-Aggregation algorithm. First, Micro-Aggregation in HUT
does not impact the DP protection since K-Means Clustering constrains only
small-valued data; and second, Micro-Aggregation in HUT does not change the
data, and therefore no negative impacts are incurred on DP protection.

Fig. 2: A demo of the two different data set settings. The first situation is
querying speed values from the Image-Speed Pair which simulates Simple Query,
while the second special situation is querying the counting numbers from Speed
Value Counts, which simulates Counting Query.

##  3 Experimental Setup

$\bullet$ Selected Dataset: We use BROOK dataset [14]. BROOK dataset is a
multi-modal dataset with facial videos and 11 dimensions of
biosignal/contextual data from 34 drivers’ driving experiences. In our
experiments, we specifically handcraft a dataset that contains drivers’ facial
images and the corresponding vehicle speeds, as a proof-of-concept. Figure 2
demonstrates an example of our handcrafted dataset.

$\bullet$ Configurations of Proposed Methods: There are three key parameters
for our method, the scale of DP $\epsilon$, the threshold $p$ and the number
of clusters $k$ for Micro-Aggregation. We testify our design in various
parameter settings, including $\epsilon\in\\{0.008,0.01,0.02,0.05\\}$,
$k\in\\{5,8,10,15\\}$ and $p\in\\{30$%$,35$%$,40$%$\\}$, in which
$k\in\\{5,8,10\\}$ are used for simple queries and $k\in\\{5,10,15\\}$ are
used for counting queries to satisfy different data set configurations. We
conduct 20 experiments for each parameter setting and take the mean as the
final value of information loss.

$\bullet$ Query Settings: We examine the data utility of the DP-protected
dataset with two representative types of queries. The first query type is a
simple query with facial images as input, and a query for the speed data in
our handcrafted dataset. These queries simulate the unit-length querying
process for intelligent in-vehicle systems to continuously retrieve data for
decision-making in fine granularity. The second type of query is the counting
query, which simulates the statistical analysis function of in-vehicle
systems. The counting query is similar to query data from a contingency table
(the bottom image in Figure 2). Different from the first general situation,
this one already has a fixed sensitivity of 1, making it more challenging for
HUT because sensitivity reduction is limited.

$\bullet$ Evaluation Metric: We use the Mean Square Error (MSE) as the
evaluation metric of data utility and information loss. We denote the query
response on the DP-protected data as $d^{\prime}$, and the query response on
the raw data as $d$. MSE is calculated by equation
$\mathrm{MSE}=\frac{1}{n}\sum_{i=1}^{n}\left(d^{\prime}_{i}-{d}_{i}\right)^{2}$.

##  4 Experimental Results

We compare HUT with the state-of-the-art DP algorithms to examine its
effectiveness of data utility. Since the average information loss of method ➋
and ➌ in Figure 1 are highly similar under all experiment conditions, we
select the overall best performer ➌ to compare its information loss with the
information loss of HUT.

The results of the simple query and counting query in different $\epsilon$
values with fixed k clusters and threshold portions are presented in Figure 3
and Figure 4. Both figures prove the effectiveness of Micro-aggregation in
increasing SNR for small-value data and mitigating the negative effects of DP
protections, in a common setting (simple query) and challenging setting
(counting query). For simple query, our method achieves a maximum of 98% (when
$\epsilon$=0.008) and a minimum of 79% (when $\epsilon$=0.05) reduction
information loss reduction compared to the state-of-the-art, with k=10 and 30%
threshold. For the counting query, the max information loss reduction
percentage is 65% (when $\epsilon$=0.008) and the minimum percentage is 27%
(when $\epsilon$=0.05). We note that when the noise magnitude grows (i.e.
$\epsilon$ values decreases), the performance of our algorithm degrades in a
simple query, but improves in counting query. The opposite trend in two query
settings proves two things: 1) Applying DP on raw data could significantly
harm normal data utility, as simple query, a primary indexing operation could
cause relatively severe information loss; 2) Aggregation, as one kind of
counting operations, could effectively avoid over-protection in extreme
conditions, but is less effective when data is not over-protected.

Fig. 3: Average information loss of different epsilon values for simple
queries (when k = 10, threshold = 35$\%$).  Fig. 4: Average information loss
of different epsilon values for counting queries (when k = 10, threshold =
35$\%$).  Fig. 5: Average information loss of different k values for simple
queries (when epsilon = 0.02, threshold = 35$\%$).  Fig. 6: Average
information loss of different k values for counting queries (when epsilon =
0.02, threshold = 35$\%$).

Figure 5 and Figure 6 presents the information loss of simple query and
counting query in different K settings, when $\epsilon$ and threshold are
fixed at 0.02,35%. We obtain the key observations that an effective choice of
the number of clusters for K-mean clustering falls in the range between 5 and
15, which have the optimal information loss reduction. In addition, our
empirical study shows that any reasonable threshold values would incur
negligible perturbations on the information loss, as different threshold
values cause a maximum 1% variance for simple query and 6% for counting query.

Table 1 reports the best percentage of average information loss compared to
the state-of-the-art mechanism in both simple query and counting query for
each $\epsilon$ value when the number of clusters k and threshold for
aggregation varies for each selected $\epsilon$. We believe exploiting the
mitigation of DP’s negative effects in different $\epsilon$ settings is
essential, as a small scale of $\epsilon$ changes could cause relatively large
fluctuations in the information loss.

Table 1: Best information loss reduction on the state-of-the-art mechanism in different $\epsilon$ settings | Epsilon | K | Threshold | % of reduction
---|---|---|---|---
Simple Query | 0.008 | 10 | 30% | 78.51%
0.01 | 8 | 30% | 80.74%
0.02 | 8 | 30% | 88.75%
0.05 | 5 | 30% | 95.69%
Counting Query | 0.008 | 10 | 35% | 64.13%
0.01 | 5 | 30% | 71.71%
0.02 | 5 | 35% | 45.22%
0.05 | 10 | 40% | 50.46%

##  5 Conclusions and Future Works

We present HUT, a new algorithm to enable High UTility for DP-enabled
protection in the context of IoV. Our key insight is to take advantage of the
fact that, the interactions between centralized servers and edge vehicles are
usually frequent and fine-grained, which are combined into a series of
unbalanced batches. We evaluate the effectiveness of HUT against the state-of-
the-art DP protection mechanism, and the results show that HUT can provide
much lower information loss by 95.69% and simultaneously enable strong
mathematically-guaranteed protection of sensitive data.

Our future works aim to explore the scalability and robustness of HUT. For
scalability, we expect to examine whether HUT can be applied for data with
broad range or unevenly distributed datasets; as for robustness, HUT still
suffers from the instability of the K-mean clustering method, which can
produce query results with relatively large variance. Therefore, we can expect
a trade-off between different clustering algorithms, in terms of the
robustness of HUT.

## References

  * [1] Qian Mei, Hu Xiong, Yanan Zhao, and Kuo-Hui Yeh,  “Toward blockchain-enabled iov with edge computing: Efficient and privacy-preserving vehicular communication and dynamic updating,”  in 2021 IEEE Conference on Dependable and Secure Computing (DSC), 2021, pp. 1–8.
  * [2] Kai Fan, Wei Jiang, Qi Luo, Hui Li, and Yintang Yang,  “Cloud-based rfid mutual authentication scheme for efficient privacy preserving in iov,”  Journal of the Franklin Institute, vol. 358, no. 1, pp. 193–209, 2021.
  * [3] Zeyu Xiong, Jiahao Wang, Wangkai Jin, Junyu Liu, Yicun Duan, Zilin Song, and Xiangjun Peng,  “Face2Statistics: User-Friendly, Low-Cost and Effective Alternative to In-Vehicle Sensors/Monitors for Drivers,”  in International Conference on Human-Computer Interaction, 2022\.
  * [4] Zhentao Huang, Rongze Li, Wangkai Jin, Zilin Song, Yu Zhang, Xiangjun Peng, and Xu Sun,  “Face2Multi-modal: In-vehicle Multi-modal Predictors via Facial Expressions,”  in Adjunct Proceedings of the 12th International Conference on Automotive User Interfaces and Interactive Vehicular Applications, AutomotiveUI 2020, Virtual Event, Washington, DC, USA, September 21-22, 2020. 2020, pp. 30–33, ACM.
  * [5] Jiahao Wang, Zeyu Xiong, Yicun Duan, Junyu Liu, Zilin Song, and Xiangjun Peng,  “The Importance Distribution of Drivers’ Facial Expressions Varies over Time!,”  in AutomotiveUI ’21: 13th International Conference on Automotive User Interfaces and Interactive Vehicular Applications, Leeds, United Kingdom, September 9-14, 2021 - Adjunct Proceedings. 2021, pp. 148–151, ACM.
  * [6] Yu Zhang, Wangkai Jin, Zeyu Xiong, Zhihao Li, Yuyang Liu, and Xiangjun Peng,  “Demystifying Interactions Between Driving Behaviors and Styles Through Self-clustering Algorithms,”  in International Conference on Human-Computer Interaction, Heidi Krömker, Ed., 2021.
  * [7] Cynthia Dwork,  “Differential privacy,”  in Automata, Languages and Programming, Michele Bugliesi, Bart Preneel, Vladimiro Sassone, and Ingo Wegener, Eds., Berlin, Heidelberg, 2006, pp. 1–12, Springer Berlin Heidelberg.
  * [8] Fiodar Kazhamiaka, Matei Zaharia, and Peter Bailis,  “Challenges and opportunities for autonomous vehicle query systems,”  in 11th Conference on Innovative Data Systems Research, CIDR 2021, Virtual Event, January 11-15, 2021, Online Proceedings. 2021, www.cidrdb.org.
  * [9] Yicun Duan, Junyu Liu, Wangkai Jin, and Xiangjun Peng,  “Characterizing Differentially-Private Techniques in the Era of Internet-of-Vehicles,” Technical Report-Feb-03 at User-Centric Computing Group, University of Nottingham Ningbo China, 2022.
  * [10] Michael Hay, Vibhor Rastogi, Gerome Miklau, and Dan Suciu,  “Boosting the accuracy of differentially private histograms through consistency,”  Proc. VLDB Endow., vol. 3, no. 1, pp. 1021–1032, 2010.
  * [11] Jordi Soria-Comas, Josep Domingo-Ferrer, David Sánchez, and Sergio Martínez,  “Enhancing data utility in differential privacy via microaggregation-based k -anonymity,”  VLDB J., vol. 23, no. 5, pp. 771–794, 2014.
  * [12] Bo-Chen Tai, Szu-Chuang Li, and Yennun Huang,  “K-aggregation: Improving accuracy for differential privacy synthetic dataset by utilizing k-anonymity algorithm,”  in 31st IEEE International Conference on Advanced Information Networking and Applications, AINA 2017, Taipei, Taiwan, March 27-29, 2017, Leonard Barolli, Makoto Takizawa, Tomoya Enokido, Hui-Huang Hsu, and Chi-Yi Lin, Eds. 2017, pp. 772–779, IEEE Computer Society.
  * [13] S. Gopal Krishna Patro and Kishore Kumar Sahu,  “Normalization: A preprocessing stage,” 2015.
  * [14] Xiangjun Peng, Zhentao Huang, and Xu Sun,  “Building BROOK: A multi-modal and facial video database for human-vehicle interaction research,”  CoRR, vol. abs/2005.08637, 2020.