11institutetext: Faculty of Informatics, Masaryk University, Botanická 68a,
602 00 Brno, Czech Republic 11email: {olha, 492606, xslanin,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Learned Indexing in Proteins: Extended Work on Substituting Complex Distance
Calculations with Embedding and Clustering Techniques

Jaroslav Olha  0000-0003-1824-468X Terézia Slanináková  0000-0003-0502-1145
Martin Gendiar  Matej Antol  0000-0002-1380-5647 Vlastislav Dohnal
0000-0001-7768-7435

###### Abstract

Despite the constant evolution of similarity searching research, it continues
to face the same challenges stemming from the complexity of the data, such as
the curse of dimensionality and computationally expensive distance functions.
Various machine learning techniques have proven capable of replacing elaborate
mathematical models with combinations of simple linear functions, often
gaining speed and simplicity at the cost of formal guarantees of accuracy and
correctness of querying.

The authors explore the potential of this research trend by presenting a
lightweight solution for the complex problem of 3D protein structure search.
The solution consists of three steps – (i) transformation of 3D protein
structural information into very compact vectors, (ii) use of a probabilistic
model to group these vectors and respond to queries by returning a given
number of similar objects, and (iii) a final filtering step which applies
basic vector distance functions to refine the result.

††The publication of this paper and the follow-up research was supported by
the ERDF ”CyberSecurity, CyberCrime and Critical Information Infrastructures
Center of Excellence” (No.CZ.02.1.01/0.0/0.0/16_019/0000822). Part of this
work was carried out with the support of ELIXIR CZ Research Infrastructure (ID
LM2018131, MEYS CR). Computational resources were supplied by the project
”e-Infrastruktura CZ” (e-INFRA CZ LM2018140 ) supported by the Ministry of
Education, Youth and Sports of the Czech Republic.

###### Keywords:

protein database, embedding non-vector data, learned metric index, similarity
search, machine learning for indexing

##  1 Introduction

The methods and approaches developed by the similarity searching community are
used by a wide range of scientific fields, both within computer science and
beyond. While some applications require provable accuracy guarantees,
verifiable algorithms, or support for complex similarity functions, many of
the similarity searching problems emerging in various areas of research do not
have such strict formal constraints.

Often, the representations of data or their similarity functions are
inherently unreliable, so the idea of a perfect similarity search is
unattainable regardless of the quality of the search method. Other times, the
problems are too computationally complex, and instead of a perfect answer that
would take too long to compute, a best-effort answer in a reasonable amount of
time is preferable (perhaps to be refined later using more costly methods).
Since the relevance of research on similarity searching is determined by its
ability to solve current problems, it should strive to employ the widest
possible array of problem-solving tools.

There are multiple similarity search indexes, falling under the umbrella of
approximate searching, capable of adjusting to such use cases by lowering
their accuracy thresholds or returning partial results [22]. However, in
recent years, an entirely new approach has begun to gain traction – along with
most other areas of computer science, the area of data retrieval has started
to incorporate various machine learning approaches. These approaches have the
potential to not only expand the available toolbox of problem-solving
instruments, but also to offer a new way of thinking about the problems
themselves. Notably, in 2018, Kraska et al. [11] suggested that all
conventional index structures could be viewed as models of data distributions,
implying that machine and deep learning models could be used in their place.
Even though the idea was originally proposed and tested on structured data,
this reframing of the problem has already inspired similar work in the realm
of unstructured datasets [1][7][21].

To investigate the potential of these approaches further, we have chosen to
examine the problem of 3D protein structure similarity search [16]. This is an
important open problem in biochemical and medical research, which can be
viewed as an instance of similarity searching in non-vector datasets, because
similarity between a pair of protein structures is usually calculated using a
series of non-trivial, computationally expensive operations. Additionally, the
amount of 3D protein structure data is currently exploding due to a recent
breakthrough in the field [9], and the demand for versatile similarity
searching approaches is likely to grow in the near future.

In this paper, we demonstrate that even a relatively complex interdisciplinary
problem such as 3D protein structure retrieval can be tackled with fast and
lightweight solutions. We present a simple pipeline where protein structures
are first transformed into short vectors and used to train multiple
partitioning and classification models – these are linked together to form a
learned index structure. The index then answers queries by returning several
candidate leaf nodes, and filtering the objects stored therein using basic
vector (similarity) functions111The entire functionality is publicly available
as a search engine prototype at https://disa.fi.muni.cz/demos/lmi-app-
protein/.

This approach, while based on probabilities rather than mathematical
guarantees, provides a reasonable quality of results at a fraction of the
computational costs required by previous methods. In addition, its modularity
allows us to change algorithms or their parameters for various trade-offs
between complexity and accuracy, depending on the particular use case.

##  2 Related work

Learned indexing was first introduced in [11] with the core idea of learning a
cumulative distribution function (CDF) to map a key to a record position. This
proposition challenged the long-standing paradigm of building indexes solely
with classic data structures such as B-trees and allowed for reduction in
searching costs compared to the traditional methods. To allow for indexing of
large data collections, the authors introduced Recursive model index (RMI) – a
hierarchical tree index structure of simple interconnected machine learning
models, each learning the mapping on a subset of the dataset. RMI is, however,
limited to sorted datasets of structured data, and cannot accommodate multi-
dimensional data.

The generalization of the learned indexing concept to spatial and low-
dimensional data was explored primarily by the Z-order model [25], which makes
use of the space-filling curve encoding to produce 1D representation of the
original data, and ML-index [6] which achieves the same with the use of
iDistance [8]. RSMI [17] introduced a recursive partitioning strategy,
building on the Z-order model. Furthermore, LISA [13] and Flood [15] both
partition the data into grid cells based on data distribution, improving the
range and kNN performance of prior approaches. These approaches were recently
directly compared to (and surpassed by) LIMS [21], which generalizes to metric
spaces and establishes a new state-of-the-art performance on datasets of up to
64 dimensions.

Indexing solutions for approximate, rather than precise queries were explored
by Chiu et al. [5] introducing a probability-based ranking model able to learn
the neighborhood relationships within data, which can be integrated with
existing product quantization (PQ) methods. This method was tested on
1-billion datasets with approximately 100 dimensions and has been shown to
boost the performance of traditional PQ methods.

Following the architectural design of RMI, we proposed the Learned metric
index (LMI) [1], which can use a series of arbitrary machine learning models
to solve the classification problem by learning a pre-defined partitioning
scheme. This was later extended to a fully unsupervised (data-driven) version
introduced in [19], which is utilized in this work.

Finally, the ideas of creating an indexing solution with machine learning have
found their use in many different domains, e.g., in trajectory similarity
search [18] or information retrieval, using a single transformer model [20].

### Protein representation

To enable computational approaches to the problem of uncovering functional
properties of proteins, a great amount of research attention has been directed
to creating representative (numerical) embeddings of protein structures. There
are two distinct categories of embeddings based on the input data – those that
operate with sequences and those working with 3D structures. Although they
serve a similar purpose, these two categories are completely distinct in terms
of the technical approaches they employ. Specifically, sequence embeddings use
techniques such as hidden Markov models [10] or various natural language
processing techniques [2] to derive meaning from protein sequences, treating
them as encoded sequences of characters – this approach is not applicable to
our research.

Embeddings representing protein 3D structures are generally less elaborate,
since the information content is more robust to begin with. The most common
encoding is a protein distance map, which produces a symmetric 2D matrix of
pairwise distances between either atoms, groups of atoms or amino acid
residues. This distance map can be transformed into a protein contact map,
which is a binary image where each pixel indicates whether two residues are
within a certain distance from one another or not. Contact maps have been used
in conjunction with machine learning techniques for prediction of protein
structure [4]; another studied problem is the reconstruction of 3D structure
based on information in contact maps [24]. While these techniques are related
to our own approach, we produce embeddings that are considerably more compact
and reflective of our similarity searching use case, as will be shown in the
following sections.

##  3 Data domain

We have chosen to test our approach on 3D protein structures for several
reasons. First, while protein structure data is very widely used, and the
study of this data is vital for almost every area of biochemical research, the
issue of efficient search and comparison of protein structures is still
unresolved to some extent, with many databases still relying on time-consuming
brute-force linear search [14]. This data is also publicly available in a
single database, called the Protein Data Bank (PDB), which is used by the
majority of protein researchers and widely agreed upon as the standard. Even
though this database is large enough to require efficient search methods, its
size ($\sim 500,000$ structures as of 2022) still makes it possible to
download the entire database and scan it exhaustively for ground truth answers
if necessary.

Just as importantly, it is clear that the issue of efficient search within
this data will only become more crucial and challenging in the next few years.
Firstly, the common dataset of empirically solved protein structures continues
to grow exponentially [3] – its current contents constitute a mere fraction of
all the protein structures in nature, and the complex laboratory procedures
needed to obtain these structures are being refined every year. Secondly, the
computational prediction of protein structures from their sequences has
recently seen rapid improvement with the release of AlphaFold 2 in 2021 [9] –
this has already resulted in hundreds of thousands of new reliable 3D
structures that are of great interest to researchers, with tens of millions
more to be added in the coming months [23].

Interestingly, the question of efficient protein structure search is not only
important for storage and retrieval purposes – for instance, researchers often
discern the function of unknown proteins by comparing them to other, better-
known proteins. Since the function of a protein is entirely dictated by its 3D
structure, this is equivalent to searching a database for the most similar
protein structures. However, the specific needs of this type of research can
vary – while some researchers are looking for extremely specific deviations
among a group of very similar proteins (e.g., when studying mutations or
conformational changes), others might be looking for much broader patterns of
similarity between distant protein families.

Structurally, a protein chain consists of a linear sequence of interconnected
building blocks called amino acids. Within the right biochemical environment,
this linear sequence folds in on itself to form a complex 3D structure.
Protein structures are sometimes cited as a typical example of complex
unstructured data, since they cannot be meaningfully ordered according to any
objective criteria (any search method needs to rely on pairwise similarity),
and the similarity of two protein structures often cannot be determined by a
single vector operation. Typically, protein molecule data are stored using the
three-dimensional coordinates of each of their atoms, with the protein
randomly oriented in space. In order to compare a pair of protein structures,
they first need to be properly spatially aligned in terms of translation and
rotation, and a subset of atoms must be selected for alignment. This typically
involves gradual optimization of a spatial distance metric (such as the root-
mean-square deviation of all the atom coordinates), which is a computationally
expensive process that cannot be directly mapped to a simple vector operation.

One commonly-used measure of protein similarity is the $Q_{score}$ [12], which
is calculated by dividing the number of aligned amino acids in both protein
chains by the spatial deviation of this alignment and the total number of
amino acids in both chains. Even though this measure is imperfect and not
appropriate for all use cases, it is used in several prominent applications,
including the PDB’s own search engine.

Note that two identical structures have a $Q_{score}$ of 1, and completely
different structures have a Q-score of 0: as a result, the score needs to be
inverted in order to be used as a distance metric ($d(x,y)=1-Q_{score}(x,y)$).
In the following sections, we will refer to this inverted value as
$Q_{distance}$.

While there are a few types of protein structure embeddings that are invariant
to the spatial alignment of the molecules (see the Related Work section),
these were not developed for the purpose of fast data retrieval. As a result,
they tend to be too detailed and cumbersome to serve as simple data
descriptors. By contrast, the embedding presented in this work has been
designed specifically to contain the optimal amount of information for
efficient and accurate similarity searching, as will be shown in the following
sections.

##  4 Fast searching in proteins

Figure 1: A diagram of the proposed solution.

We present a pipeline (visualized in Figure 1) consisting of three separate
components: (i) a simple embedding technique for protein data in the PDB
format, (ii) the use of a machine-learning-based index – Learned Metric Index
(LMI) – to locate a candidate set of similar protein structures, and (iii)
fast filtering to produce the final query answer.

The embedding we propose divides the protein sequence into $N$ consecutive
sections – the positions of the atoms within each section are averaged, and
the section is subsequently treated as a single point in space. We then
calculate distances between each pair of these sections, creating an incidence
matrix. In this matrix, we prune all the values exceeding a cutoff, and
normalize the rest. The matrix is symmetrical and all the diagonal values are
0. The half of this matrix (omitting the main diagonal) is then reduced into a
single row in a $M\times(\frac{N^{2}-N}{2})$ matrix, where $M$ is the number
of proteins in the entire dataset (see Figure 1).

This produces a very compact embedding for all the proteins, and the entire
dataset can be represented by a file that is up to two orders of magnitude
smaller than the original database – see Table 1.

Embedding Size | File Size | Index Build Time (256-64) | Index Build Time (128-128)
---|---|---|---
$5\times 5$ |  $16$ MB | 246s | 184s
$10\times 10$ |  $51$ MB | 350s | 270s
$30\times 30$ |  $456$ MB | 927s | 655s
$50\times 50$ |  $1275$ MB | 2391s | 1814s
Table 1: File size of the protein dataset (518,576 protein chains) stored
using protein embeddings, and build times of two different LMI architectures.
Note that the size of the original database is 8.2 GB.

To reduce the search space to a small number of candidate protein structures,
we used the Learned Metric Index (LMI), a tree index structure where each
internal node is a learned model trained on a sub-section of data assigned to
it by its parent node [1]. Specifically, we used the data-driven version of
LMI, where the partitioning is determined in an unsupervised manner. We
explored different architectural setups – both in terms of the number of nodes
at each level (index breadth), as well as the number of levels (index depth).
As the learned models, we explored K-Means, Gaussian Mixture Models, and
K-Means in combination with Logistic regression (see [19] for details
regarding the model setups). For the sake of compactness, in the experimental
evaluation we only present the results achieved with the best-performing setup
– a two-level LMI structure with arity of 256 on level 1 and 64 on level 2
(i.e., 256 root descendants, each of them with 64 child nodes), with K-Means
chosen as the partitioning algorithm. After LMI is built, we search within it
using a query protein structure and return target candidate sets; the size of
the candidate sets is determined by a pre-selected stop condition (for
instance, a stop condition of 1% of the dataset corresponds to $\sim 5,000$
candidate answers per query).

In the final step, we filter the candidate set according to a particular
distance function. In our experiments, we have examined filtering based on the
Euclidean distance as well as the cosine distance of the vector embeddings,
but the filtering step could theoretically be performed using any distance
metric, or even the original $Q_{score}$ similarity of the full protein
structures. The filtering step returns a subset of the candidate set based on
the specified criteria (i.e., kNN or range).

##  5 Experimental evaluation

We evaluated our approach using range queries, with 512 randomly chosen
protein chains from the dataset used as query objects. In order to compare our
results against the ground truth, we needed to know the $Q_{distances}$ (based
on $Q_{score}$) between the 512 protein chains and all the other chains in the
database – these distances were kindly provided by the researchers behind
[14], where the same 512 objects were used as the pivots for their search
engine. The objects were chosen uniformly randomly with respect to protein
chain length, which ensures that even very long proteins are represented among
our queries (despite constituting a relatively small portion of the dataset).

Figure 2: Evaluation of range queries after LMI search and before filtering,
using the K-Means method and a 256-64 LMI architecture: (left) Range=0.1,
(middle) Range=0.3, (right) Range=0.5.

We expected the performance of our method to deteriorate as the range of the
queries expands, since a wider search range would require the method to
correctly identify more objects which are less similar. To examine this
effect, we have chosen three representative query ranges of 0.1, 0.3 and 0.5 –
in a real use case, the range would be chosen by a domain expert based on the
particular use case. As a rule of thumb, a range of 0.1 represents a high
degree of similarity, while a range of 0.5 represents low (but still
biologically significant) similarity; the biological relevance of answers
drops sharply beyond this range [14].

First, we evaluated the performance of the LMI, before the filtering step. The
recall shown in Figure 2 pertains to the entire candidate set of objects (i.e.
how much of the ground truth answer is contained in the 1%/5%/10% of the
dataset returned by the LMI for further filtering)222Precision is not
evaluated in this step – at this point in the pipeline, it is very low and not
particularly relevant..

This figure presents us with two important pieces of information – firstly, it
is clear that LMI can reach very high recall even when trained on the smaller
10x10 embedding – this makes the embedding a natural choice for further
evaluation, since it is efficient while significantly reducing the memory and
CPU costs of training compared with the larger embeddings. It can also be seen
that, especially in the lower query ranges which are of most interest to us,
the 1% stop condition represents a sensible trade-off between recall and
search time, returning relatively few candidate objects ($\sim 5,000$) while
minimizing the amount of false negatives.

Figure 3: Distribution of objects in the LMI leaf nodes. (A completely
balanced structure would hold $\sim 30$ objects in each bucket).

Figure 4: Correlation between $Q_{distance}$ and the Euclidean distance used
for filtering: Blue points represent proteins returned by the LMI; orange
points represent proteins which are present in the ground truth answer but
were not returned by the LMI.

In addition to recall, we also need to investigate the distribution of objects
in the index, to ensure that the occupancy of leaf nodes (i.e. data buckets)
is not overly imbalanced – an extremely imbalanced index would achieve high
recall, but it would also be more likely to return overly large candidate
sets, which would be detrimental to the filtering step. The size distribution
of the buckets is shown in Figure 4, and it confirms that the distribution is
not overly skewed towards large buckets, even when the embeddings are quite
small, as is the case with N=10. However, the embeddings can only be condensed
up to a point – the 5x5 embedding (which transforms each protein structure
into a vector of only 10 values) causes a large portion of the objects to
concentrate in a single bucket, as the LMI can no longer distinguish among
groups of objects.

During the filtering step, recall naturally decreases over time (since the
method occasionally filters out relevant answers), while precision should
improve as the portion of relevant objects in the candidate set increases.
This effect will be different depending on the distance metric used for
filtering, as well as the dataset – we show the effects for two distance
metrics (Euclidean distance and cosine distance) in Figure 5. By analyzing the
correlation between $Q_{distance}$ and each of the two distance metrics, we
have determined that the Euclidean distance is the better filtering function
for this dataset.333Note that the actual $Q_{distance}$ function scales
differently from the Euclidean distance used to filter answers (see Figure 4)
– this requires simple re-scaling (e.g., to find range=0.5 queries, we set the
Euclidean distance cutoff at 0.75).

Figure 5: Effects of filtering on the recall and precision of the candidate
set of objects (relative to the ground truth answer).

Table 2 shows the final results of the range queries with the best-performing
configuration of parameters: embedding size N=10, the K-means clustering
model, 256-64 LMI architecture, and filtering after the 1% stop condition
using the Euclidean distance metric. The results, especially in the lower
query ranges, are very encouraging, although the filtering stage seems to
introduce a surprisingly large amount of false negatives by filtering out
parts of the correct answer. It is likely that the filtering metric we have
chosen was slightly too naïve, and the filtering step could have benefited
from a different distance function, or at least a different weighting of the
vectors before calculating their Euclidean distance. In the future, this
presents a natural point of focus to improve our results even further.

| Range 0.1 | Range 0.3 | Range 0.5
---|---|---|---
| Mean # of objects: 83 | Mean # of objects: 236 | Mean # of objects: 519
LMI Recall | 0.973 (1.000) | 0.895 (0.999) | 0.755 (0.867)
Recall after filtering | 0.742 (0.878) | 0.649 (0.711) | 0.530 (0.637)
F1 after filtering | 0.712 (0.855) | 0.669 (0.766) | 0.592 (0.673)
Table 2: Overall evaluation of protein range queries: the average values, as
well as the median values (in parentheses) are shown.

Figure 4 further illustrates the strengths and weaknesses of our approach by
showing the relationship between the $Q_{distance}$ metric and the Euclidean
distances of vector embeddings. There is clear correlation between these two
metrics, making a strong case for the simpler one, and LMI is much more
successful at finding objects that are more similar to the query (left side of
the graph) than it is at finding less similar objects (right side of the
graph). It should be restated that even though we use the $Q_{distance}$
metric as the ”ground truth”, it is merely a subjective similarity metric, and
the relevance of results that are close to the edge of the query range is
debatable (i.e., an object with a $Q_{distance}$ of 0.499 is not necessarily
more relevant than an object with a $Q_{distance}$ of 0.501).

Distance computations on long proteins pose a considerable challenge for
similarity searching methods, typically requiring a disproportionate amount of
computing time [14]. This does not apply to our embedding approach, which
transforms all the proteins into fixed-length vectors using the same
algorithm, regardless of their length. It would therefore intuitively follow
that such an approach should achieve lower recall on the longer protein
queries, since more information is lost in the embedding. However, in
practice, that is not the case, as can be seen in Figure 6.

Figure 6: Distribution of recall for different protein chain lengths (left to
right: the shortest 10% of chains, each quartile from shortest to longest, and
the longest 10% of chains).

This is probably due to the simple fact that the distribution of chain lengths
in protein databases is not uniform – in fact, the long chains are much less
numerous. As a result, losing a lot of information about the long chains is
not a problem, since they are still relatively easy to find. This gives the
fixed-length embedding approach a significant performance advantage – in a
dataset where there are relatively few long protein chains, it would be a
waste of resources to calculate and store an excessive amount of data about
them, the way many other methods do.

As has been mentioned in the previous section, in addition to variable chain
lengths, the proteins also have a variable number of neighbors in any given
range. This can, again, call into question the reliability of the results
based on recall, since recall is a relative measure, and the size of the error
will be different based on the size of the actual answer. For instance, if all
range(0.1) query answers only consisted of a single easy-to-find object,
evaluating recall on its own would give us a biased idea of the searching
efficiency. Figure 7 shows how many neighbors each protein structure has
according to the ground truth, and how many of these neighbors have been found
using our approach – while there are significantly more proteins with fewer
than 10 neighbors in the lower query ranges, these are not the majority, and
the errors are distributed evenly relative to query answer size.

Figure 7: The number of results returned by the LMI compared to the ground
truth. The X-axis shows the absolute number of objects in a given query answer
(according to the ground truth), while the Y-axis shows the portion of the
objects that were found by LMI (i.e., the recall). Each blue point in the
scatter plot represents one query, with the larger and more warmly-colored
points representing multiple queries with identical graph coordinates.

Finally, to provide broader context for the pipeline’s performance, we have
evaluated it against a more conventional approach. While there is no analogous
method for searching 3D protein structures using range queries, there is a
similar, recently-published method for nearest neighbor retrieval in the same
database of 3D protein structures [14]. This method uses a three-stage search
engine which compares bit-strings in the Hamming space (”sketches”) to
approximate the distance of protein chains.

To allow for a relatively fair comparison of these two methods, we modified
our search parameters to more closely match the search results presented in
[14]. Specifically, since their similarity queries were mainly 30NN queries
limited by the range 0.5, we have performed range queries in the same range
(even though our method’s performance is substandard in this range), and in
the filtering step, we filtered out all objects beyond the 30 best-ranking
answers. Since the sketch-based method was originally compared with the linear
search of the PDB database, we present this benchmark as well – naturally,
this is always the slowest method by far, but it requires no index and always
finds the exact answer.

All of these results can be found in Table 3 – while our method clearly does
not match the high accuracy of the sketch-based method in this experimental
setup, it is faster by at least an order of magnitude, occasionally exceeding
4 orders of magnitude since it does not suffer from an extreme ”tail” of
worst-case search times caused by evaluation of long proteins.

| LMI + Filtering 444Before the filtering step, the candidate set returned by LMI has a median recall of 1.0 and an average recall of 0.87; however, since the candidate set is much larger than the final answer ($\sim 5,000$ objects), the accuracy is insignificant and has thus been omitted from the table. | Sketch-based method | PDB Engine
---|---|---|---
Accuracy (median) | $0.660$ | $1.0$ |  $1.0$
Accuracy (mean) | $0.626$ | $0.937$ |  $1.0$
Time (median) |  $0.094$ s |  $2.5$ s |  $183$ s
Time (max) |  $0.145$ s |  $6109$ s |  $14321$ s
Index size555The size of the internal structure of the index, excluding raw protein data. |  $87$ MB |  $178$ MB 666495,085 * (320b (sketches) + 1024b (sketches) + 4*6*64b (PPP-codes)) + 512 * 16kB (pivots) | N/A
Table 3: The accuracy, search times, and memory requirements of 30NN protein
search queries with a maximum distance radius of 0.5.

##  6 Summary and Conclusions

In an effort to investigate the potential of new data retrieval techniques in
the field of similarity searching, we have developed and evaluated a novel
approach to the problem of protein structure search, resulting in a short
pipeline consisting of a concise vector embedding, learned indexing, and
distance-based answer filtering. By successfully applying this approach on a
well-established database of 3D protein structures, we have shown that even in
a domain that may, at first, seem poorly suited to simple vector-based
transformations, a surprising amount of information can be discerned by
learned models.

One advantage of our modular approach is that every part of the pipeline can
be evaluated separately, allowing experts to identify the weakest spots and
alter them based on the current use case and dataset. The experiments
presented in this paper serve as a good example – after evaluating each part
of our own pipeline, it is clear that we have chosen an overly simplistic
filtering method for our data. In the future, we plan to investigate more
sophisticated options for vector-based filtering, as well as a completely
different approach to reducing the size of query answers.

While it is difficult to compare our work with the state of the art (since
there are no direct analogues to our method in the chosen data domain), we
have made an effort to modify our method for the fairest possible comparison
with a recent, more conventional similarity searching approach in the same
domain. In this comparison, our solution, although coming up short in terms of
accuracy, is consistently faster by multiple orders of magnitude, and
maintains much lower memory requirements.

Our work aims to make a case for less conventional solutions to similarity
searching problems – ones that rely on learned approaches, rather than
traditional indexing models, to discern the natural distribution of data. This
is by no means an argument for universal adoption of machine-learning-based
techniques in all similarity searching applications. The approach in this
paper is problem-specific, and has several potential drawbacks that need to be
weighed against its strengths. Nevertheless, we provide some insight as to how
the trend of machine-learned pattern recognition, which has already caused
minor and major revolutions in most fields of computer science, could be
applied in practical solutions to current similarity searching problems.

Based on these results, we remain convinced that learned indexing approaches
(such as the Learned Metric Index used in this work) will play an integral
part in shaping the future of similarity searching.

## References

  * [1] Antol, M., Ol’ha, J., Slanináková, T., Dohnal, V.: Learned metric index — proposition of learned indexing for unstructured data. Information Systems 100 (2021)
  * [2] Asgari, E., Mofrad, M.R.: Continuous distributed representation of biological sequences for deep proteomics and genomics. PloS one 10(11), e0141287 (2015)
  * [3] Burley, S.K., Bhikadiya, C., Bi, C., Bittrich, S., Chen, L., Crichlow, G.V., Christie, C.H., Dalenberg, K., Di Costanzo, L., Duarte, J.M., et al.: Rcsb protein data bank: powerful new tools for exploring 3d structures of biological macromolecules for basic and applied research and education in fundamental biology, biomedicine, biotechnology, bioengineering and energy sciences. Nucleic acids research 49(D1), D437–D451 (2021)
  * [4] Cheng, J., Baldi, P.: Improved residue contact prediction using support vector machines and a large feature set. BMC bioinformatics 8(1), 1–9 (2007)
  * [5] Chiu, C.Y., Prayoonwong, A., Liao, Y.C.: Learning to index for nearest neighbor search. IEEE transactions on pattern analysis and machine intelligence 42(8), 1942–1956 (2019)
  * [6] Davitkova, A., Milchevski, E., Michel, S.: The ML-index: A multidimensional, learned index for point, range, and nearest-neighbor queries. In: EDBT. pp. 407–410 (2020)
  * [7] Hünemörder, M., Kröger, P., Renz, M.: Towards a learned index structure for approximate nearest neighbor search query processing. In: International Conference on Similarity Search and Applications. Springer (2021), to appear
  * [8] Jagadish, H.V., Ooi, B.C., Tan, K.L., Yu, C., Zhang, R.: idistance: An adaptive B+-tree based indexing method for nearest neighbor search. ACM Transactions on Database Systems (TODS) 30(2), 364–397 (2005)
  * [9] Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., Tunyasuvunakool, K., Bates, R., Žídek, A., Potapenko, A., et al.: Highly accurate protein structure prediction with alphafold. Nature 596(7873), 583–589 (2021)
  * [10] Koski, T.: Hidden Markov models for bioinformatics, vol. 2. Springer Science & Business Media (2001)
  * [11] Kraska, T., Beutel, A., Chi, E.H., Dean, J., Polyzotis, N.: The case for learned index structures. In: Proceedings of the 2018 International Conference on Management of Data. p. 489–504. SIGMOD ’18, Association for Computing Machinery (2018)
  * [12] Krissinel, E., Henrick, K.: Secondary-structure matching (ssm), a new tool for fast protein structure alignment in three dimensions. Acta Crystallographica Section D: Biological Crystallography 60(12), 2256–2268 (2004). https://doi.org/10.1107/S0907444904026460
  * [13] Li, P., Lu, H., Zheng, Q., Yang, L., Pan, G.: Lisa: A learned index structure for spatial data. In: Proceedings of the 2020 ACM SIGMOD international conference on management of data. pp. 2119–2133 (2020)
  * [14] Mic, V., Raček, T., Křenek, A., Zezula, P.: Similarity search for an extreme application: Experience and implementation. In: International Conference on Similarity Search and Applications. pp. 265–279. Springer (2021)
  * [15] Nathan, V., Ding, J., Alizadeh, M., Kraska, T.: Learning multi-dimensional indexes. In: Proceedings of the 2020 ACM SIGMOD International Conference on Management of Data. pp. 985–1000 (2020)
  * [16] Olha, J., Slanináková, T., Gendiar, M., Antol, M., Dohnal, V.: Learned indexing in proteins: Substituting complex distance calculations with embedding and clustering techniques. In: Skopal, T., Falchi, F., Lokoč, J., Sapino, M.L., Bartolini, I., Patella, M. (eds.) Similarity Search and Applications. pp. 274–282. Springer International Publishing, Cham (2022)
  * [17] Qi, J., Liu, G., Jensen, C.S., Kulik, L.: Effectively learning spatial indices. Proceedings of the VLDB Endowment 13(12), 2341–2354 (2020)
  * [18] Ramadhan, H., Kwon, J.: X-fist: Extended flood index for efficient similarity search in massive trajectory dataset. Information Sciences (2022)
  * [19] Slanináková, T., Antol, M., OǏha, J., Kaňa, V., Dohnal, V.: Data-driven learned metric index: an unsupervised approach. In: International Conference on Similarity Search and Applications. pp. 81–94. Springer (2021)
  * [20] Tay, Y., Tran, V.Q., Dehghani, M., Ni, J., Bahri, D., Mehta, H., Qin, Z., Hui, K., Zhao, Z., Gupta, J., et al.: Transformer memory as a differentiable search index. arXiv preprint arXiv:2202.06991 (2022)
  * [21] Tian, Y., Yan, T., Zhao, X., Huang, K., Zhou, X.: A learned index for exact similarity search in metric spaces. arXiv preprint arXiv:2204.10028 (2022)
  * [22] Vadicamo, L., Gennaro, C., Amato, G.: On generalizing permutation-based representations for approximate search. In: International Conference on Similarity Search and Applications. pp. 66–80. Springer (2021)
  * [23] Varadi, M., Anyango, S., Deshpande, M., Nair, S., Natassia, C., Yordanova, G., Yuan, D., Stroe, O., Wood, G., Laydon, A., et al.: Alphafold protein structure database: massively expanding the structural coverage of protein-sequence space with high-accuracy models. Nucleic acids research 50(D1), D439–D444 (2022)
  * [24] Vendruscolo, M., Kussell, E., Domany, E.: Recovery of protein structure from contact maps. Folding and Design 2(5), 295–306 (1997)
  * [25] Wang, H., Fu, X., Xu, J., Lu, H.: Learned index for spatial queries. In: 2019 20th IEEE International Conference on Mobile Data Management (MDM). pp. 569–574. IEEE (2019)


# Image Processing via Multilayer Graph Spectra

Songyang Zhang, Qinwen Deng, and Zhi Ding S. Zhang, Q. Deng, and Z. Ding are
with the Department of Electrical and Computer Engineering, University of
California at Davis, Davis, CA 95616 USA (e-mail: sydzhang, mrdeng, and
zding@ucdavis.edu).

###### Abstract

Graph signal processing (GSP) has become an important tool in image processing
because of its ability to reveal underlying data structures. Many real-life
multimedia datasets, however, exhibit heterogeneous structures across frames.
Multilayer graphs (MLG), instead of traditional single-layer graphs, provide
better representation of these datasets such as videos and hyperspectral
images. To generalize GSP to multilayer graph models and develop multilayer
analysis for image processing, this work introduces a tensor-based framework
of multilayer graph signal processing (M-GSP) and present useful M-GSP tools
for image processing. We then present guidelines for applying M-GSP in image
processing and introduce several applications, including RGB image
compression, edge detection and hyperspectral image segmentation. Successful
experimental results demonstrate the efficacy and promising futures of M-GSP
in image processing.

###### Index Terms:

Image processing, multilayer graph signal processing, convolution, sampling
theory

##  I Introduction

Signal processing over graphs has demonstrated successes in image processing
[1] owing to its ability to uncover some hidden data features. Representing
pixels (or superpixels) by nodes and their internal interactions by edges,
graphs can be formed to describe the underlying image structures. In graph
signal processing (GSP), graph Fourier space [2] and the corresponding
spectral analysis have been applied in image compression [3] and denoising
[4]. Nevertheless, single-layer graphs tend to be less efficient in processing
high-dimensional images, where each frame corresponds to a distinct geometric
structure. For example, hyperspectral images at different frequencies may
exhibit heterogeneous graph structures, and do not lend themselves to
traditional GSP tools. Similar scenarios include the color frames in RGB
images and temporal frames in videos. On the other hand, since these frames
are correlated, their inter-frame correlations should not be ignored in graph
representations. Such complex intra- and inter- frame interactions among
pixels can be represented by multilayer graphs (or multilayer networks)
instead of a single-layer graph. Shown as Fig. 1, a hyperspectral image can be
naturally modeled as a multilayer graph (MLG) with each pixel position as a
node and each frequency frame as one layer. Of strong interest are questions
of how to generalize single-layer GSP to multilayer graph signal processing
(M-GSP) and how to extract additional benefits of M-GSP in problems such as
image processing.

Graph signal processing over multilayer networks (graphs) has attracted
increasing attentions recently. In [6], a two-step graph Fourier transform
(GFT) is separately implemented spatially and temporally to compress spatial-
temporal skeleton data. Similarly, authors of [7] defined a joint time-vertex
Fourier transform (JFT) by applying GFT and DFT consecutively on time-varying
datasets. However, both the two-step GFT and JFT assume that the spatial
connections in all layers match the same underlying graph structure, thereby
limiting the use of such spectral analysis on heterogeneous data frames. In
addition to matrix-based GSP analysis, a tensor-based multi-way graph signal
processing framework (MWGSP), proposed in [8], decomposes the high-dimensional
tensor signal into individual orders and constructs factor graphs for each
order. Although MWGSP allows more complex interlayer connections and is
suitable for high-dimensional signals, it assumes the dataset to lie in a
product graph from all factor graphs, therefore still requires all intralayer
connections to exhibit homogeneous structures from the node-wise factor graph.

To capture the heterogeneous layer structures in MLG for generalized GSP, a
framework of multilayer network signal processing (M-GSP) has been proposed by
using tensor representations [9]. Defining MLG spectra via tensor
decomposition, both joint and order-wise analysis for MLG can tackle various
intralayer and complex interlayer connections. However, we note that existing
M-GSP tools still lack some definitions of basic MLG spectral operations
useful in image processing, such as sampling and convolution.

In this work, we investigate the use of M-GSP spectral analysis [9] and the
use of M-GSP tools in image processing. The contributions can be summarized as
follows:

  * •

We first introduce the MLG models and representations for images, and then
define several M-GSP tools useful for image processing.

  * •

We provide a tutorial on applying M-GSP in image processing and present
several practical cases, including image compression, edge detection, and
hyperspectral image segmentation.

We compare each proposed M-GSP image application to benchmarks from
traditional GSP and/or learning methods. Our test results demonstrate the
strength of M-GSP and the proposed MLG spectral tools in processing images.

##  II Models and Spectrum of Images in M-GSP

Figure 1: Example of Multilayer Graph Model for Hyperspectral Images: pixels
are represented by nodes and their similarities are represented by edges.

In this section, we introduce the MLG model of images and outline fundamentals
of M-GSP [9]. We refrain from using the term “network” because of its
different meanings in various field ranging from communications to deep
learning. From here onward, unless otherwise specified, we shall use the less
ambiguous term of multilayer graph instead.

To model related image frames as multilayer graph (MLG), we can naturally
define layers as frames and nodes as pixels (superpixels) in each frame.
Constructing interactions based on feature similarity or physical connections,
a multilayer graph with the same number of nodes in each layer can be
constructed to represent the images. Moreover, viewing $N$ (super)pixels
$\\{x_{1},\cdots,x_{N}\\}$ as entities, the MLG structure can be viewed as
projecting the physical entities into natural representations of $M$ frames
$\\{l_{1},\cdots,l_{M}\\}$, e.g., spectrum band frames for hyperspectral
images, color frames for RGB images, and temporal frames for time-varying
images. Such multilayer graph structure with $M$ layers and $N$ nodes in each
layer can be intuitively represented by a forth-order adjacency tensor
$\mathbf{A}=(A_{\alpha i\beta j})\in\mathbb{R}^{M\times N\times M\times N}$
with $1\leq\alpha,\beta\leq M,1\leq i,j\leq N$. Here, each entry $A_{\alpha
i\beta j}$ of the adjacency tensor $\mathbf{A}$ indicates the edge strength
between entity $j$’s projected node in layer $\beta$ and entity $i$’s
projected node in layer $\alpha$. Similar to normal graphs, Laplacian tensor
$\mathbf{L}\in\mathbb{R}^{M\times N\times M\times N}$ can be defined as the
alternative representation of MLG. Interested readers may refer to [9] for
more details of the adjacency and Laplacian tensor. For convenience, we use a
general notation $\mathbf{F}\in\mathbb{R}^{M\times N\times M\times N}$ to
represent the MLG in this work. In addition to the representation of MLG,
signals over MLG can be intuitively defined as $\mathbf{s}=({s_{\alpha
i}})\in\mathbb{R}^{M\times N}$, with $s_{\alpha i}$ as the value of pixel $i$
in frame $\alpha$.

In an undirected multilayer graph, the representing tensor
(adjacency/Laplacian) $\mathbf{F}$ is partially symmetric between orders one
and three, and between orders two and four, respectively. Then, it can be
approximated via orthogonal CANDECOMP/PARAFAC (CP) decomposition [10] as

$\displaystyle\mathbf{F}$
$\displaystyle\approx\sum_{\alpha=1}^{M}\sum_{i=1}^{N}\lambda_{\alpha
i}\cdot\mathbf{f}_{\alpha}\circ\mathbf{e}_{i}\circ\mathbf{f}_{\alpha}\circ\mathbf{e}_{i},$
(1)

where $\circ$ is the tensor outer product,
$\mathbf{f}_{\alpha}\in\mathbb{R}^{M}$ are orthonormal,
$\mathbf{e}_{i}\in\mathbb{R}^{N}$ are orthonormal and
$\tilde{\mathbf{V}}_{\alpha
i}=\mathbf{f}_{\alpha}\circ\mathbf{e}_{i}\in\mathbb{R}^{M\times N}$. Here,
$\mathbf{f}_{\alpha}$ and $\mathbf{e}_{i}$ are named as the MLG spectral bases
characterizing the features of layers and entities, respectively.

Besides MLG spectral space, the singular space is defined from HOSVD [11] as
an alternative subspace of MLG, i.e.,

$\mathbf{F}=\mathbf{S}\times_{1}\mathbf{U}^{(1)}\times_{2}\mathbf{U}^{(2)}\times_{3}\mathbf{U}^{(3)}\times_{4}\mathbf{U}^{(4)},$
(2)

where $\times_{n}$ denotes the $n$-mode product [15] and $\mathbf{U}^{(n)}$ is
a unitary $(I_{n}\times I_{n})$ matrix, with $I_{1}=I_{3}=M$ and
$I_{2}=I_{4}=N$. For an undirected multilayer graph, the adjacency tensor is
symmetric for every 2-D combination. Thus, there are two modes of singular
spectrum, i.e., $(\gamma_{\alpha},\mathbf{f}_{\alpha})$ for mode $1,3$, and
$(\sigma_{i},\mathbf{e}_{i})$ for mode $2,4$. More specifically,
$\mathbf{U}^{(1)}=\mathbf{U}^{(3)}=(\mathbf{f}_{\alpha})$ and
$\mathbf{U}^{(2)}=\mathbf{U}^{(4)}=(\mathbf{e}_{i})$. Singular tensor analysis
and spectral analysis are both efficient tools for image processing depending
on specific tasks. We shall elaborate in Section IV.

Suppose that $\mathbf{E}_{f}=[\mathbf{f}_{1},\cdots,\mathbf{f}_{M}]$ and
$\mathbf{E}_{e}=[\mathbf{e}_{1},\cdots,\mathbf{e}_{N}]$ include all the layer-
wise (frame) and entity-wise (pixel) spectral (or singular) vectors. The MLG
Fourier (or Singular) transform can be defined as
$\hat{\mathbf{s}}=\mathbf{E}_{f}^{\mathrm{T}}\mathbf{s}\mathbf{E}_{e}\in\mathbb{R}^{M\times
N}$, and the inverse transform is given as
$\mathbf{s}^{\prime}=\mathbf{E}_{f}\hat{\mathbf{s}}\mathbf{E}_{e}^{\mathrm{T}}$.
Due to page limitation, here we only introduce fundamentals of MLG spectral
and singular vectors for MLG analysis. More details, such as spectral
interpretation, frequency analysis, and filter design can be found in [9].

##  III M-GSP Tools for Image Processing

In this section, we introduce important operations for M-GSP frequency
analysis specially useful in image processing.

####  III-1 MLG Spectral Convolution

Graph-based spectral convolution is useful in applications such as feature
extraction, classification, and graph convolutional networks [12]. Graph
convolution can be implemented indirectly in the transform domain [13]. Define
graph convolution as
$\mathbf{x}\star\mathbf{y}=\mathcal{F}^{-1}(\mathcal{F}(\mathbf{x})*\mathcal{F}(\mathbf{y}))$,
where $*$ denotes Hadamard product [15], $\mathcal{F}$ is the graph Fourier
transform (GFT) and $\mathcal{F}^{-1}$ denotes inverse GFT. This formulation
exploits the property that graph convolution in vertex domain is equivalent to
spectral domain product.

In [14], a more general definition of graph convolution take the filtering
form, i.e.,
$\mathbf{x}\star\mathbf{y}=f(\mathbf{x})=(\mathbf{V}diag(\mathbf{\hat{y}})\mathbf{V}^{-1})\cdot\mathbf{x}$,
where $\mathbf{V}$ consists of eigenvectors of the representing matrix as
columns and $\hat{\mathbf{y}}$ can be interpreted as filter parameters in
graph spectral domain. Note that, these two definitions of graph convolution
have the same mathematical formulation. For convenience, we define spectral
convolution over MLG in a similar mathematical formulation as follows. In
M-GSP, the spectral convolution is defined as

$\mathbf{x}\star\mathbf{y}=\mathcal{F}_{\rm MLG}^{-1}(\mathcal{F}_{\rm
MLG}(\mathbf{x})*\mathcal{F}_{\rm MLG}(\mathbf{y})),$ (3)

where $\mathcal{F}_{\rm MLG}$ is M-GFT (or M-GST) and $\mathcal{F}_{\rm
MLG}^{-1}$ is the inverse M-GFT (or M-GST). Then, the MLG convolutional filter
with parameters $\mathbf{c}$ on the signal $\mathbf{s}$ is defined by
$f_{c}(\mathbf{s})=\mathbf{s}\star\mathbf{c}$.

###  III-A Sampling and Interpolation

####  III-A1 Sampling in Vertex Domain

Within the M-GSP framework, sampling is a process to select a subset of
signals to describe some global information. Consider the sampling of signal
$\mathbf{s}\in\mathbb{R}^{M\times N}$ over MLG into $P$ layers with index
$\\{p_{1},p_{2},\cdots,p_{P}\\}$ and $Q$ entities with index
$\\{q_{1},q_{2},\cdots,q_{Q}\\}$. Let $\times_{n}$ be the $n$-mode product
introduced in [15]. The sampling operation is defined by

$\mathbf{s}_{D}=\mathbf{s}\times_{1}\mathbf{S}_{P}\times_{2}\mathbf{S}_{Q}\in\mathbb{R}^{P\times
Q}$ (4)

where sampling operator $\mathbf{S}_{P}\in\mathbb{R}^{P\times M}$ consists of
elements

$\displaystyle[{S}_{P}]_{ij}=\left\\{\begin{array}[]{rcl}1&&j=p_{i}\\\
0&&\mbox{otherwise}\end{array},\right.$ (7)

and $\mathbf{S}_{Q}\in\mathbb{R}^{Q\times N}$ with elements calculated as

$\displaystyle[{S}_{Q}]_{ij}=\left\\{\begin{array}[]{rcl}1&&j=q_{i}\\\
0&&\mbox{otherwise}\end{array}.\right.$ (10)

The interpolation operation is defined by

$\mathbf{s}_{R}=\mathbf{s}_{D}\times_{1}\mathbf{T}_{M}\times_{2}\mathbf{T}_{N}\in\mathbb{R}^{M\times
N},$ (11)

where $\mathbf{T}_{M}\in\mathbb{R}^{M\times P}$ and
$\mathbf{T}_{N}\in\mathbb{R}^{N\times Q}$. Sampling in vertex domain can be
viewed as a two-step sampling in layers and entities. M-GSP sampling and GSP
sampling share many interesting properties. Readers may refer to [16] for more
details.

(a) Layer-Wise Sampling.

(b) Entity-Wise Sampling.

(c) Block-wise (Layer=3).

(d) Block-wise (Layer=2).

(e) Block-wise (Layer=1).

Figure 2: Results of Sampling in the Order of Transformed Value.

####  III-A2 Sampling in Spectral Domain

Beyond direct sampling in the vertex domain, it can be more convenient and
straightforward to sample signals in spectral domain. By neglecting spectral
domain zeros, we can easily approximately recover the full signal from samples
by interpolating with zeros to the sampled signal transforms before recovering
with the inverse transform. Sampling in spectral domain is intuitive for
lossless sampling and recovery. Here we shall focus on the lossy sampling of
tensor-based signals.

Generally, we can sample signals in spectral domain according to Algorithm 1.
We reshape transformed signals in terms of eigenvalues or transformed values
from top left to bottom right in descending order. When sampling a sub-block
in $\mathbf{\hat{s}}$, there are three different directions shown in Fig. 3.
Depending on the structure and features of the underlying dataset, different
methods may be suitable for different tasks. We will show examples of data
compression based on the MLG sampling theory in Section IV.

1: Input: Signals $\mathbf{s}\in\mathbb{R}^{M\times N}$;

2: Transform the signal into spectral space as
$\mathbf{\hat{s}}\in\mathbb{R}^{M\times N}$ via M-GFT or M-GST;

3: Reorder the elements in $\mathbf{\hat{s}}$ via certain rules;

4: Sample transformed signal in the left top as
$\mathbf{\hat{s}}_{D}\in\mathbb{R}^{P\times Q}$;

5: Interpolate zeros in the right bottom of $\mathbf{\hat{s}}_{D}$ as
$\mathbf{\hat{s}}_{R}\in\mathbb{R}^{M\times N}$

6: Implement the inverse MLG transform on $\mathbf{\hat{s}}_{R}$ to obtain the
recovered signal $\mathbf{s}_{R}\in\mathbb{R}^{M\times N}$.

7: Output: Sampled coefficients $\mathbf{\hat{s}}_{D}$ and recovered signal
$\mathbf{s}_{R}$.

Algorithm 1 Sampling of Tensor Signal in Spectral Domain

##  IV Application Examples

###  IV-A Image Compression

Data compression through sampling may reduce dataset size for efficient
storage and transmission. Moreover, projecting signals into a suitable
manifold or space, one may reveal a sparse representation of the data. In
M-GSP, we can also select a subset of transformed signal samples in spectral
domain to approximate the full signal as suggested in Algorithm 1. Moreover,
we can also integrate data encoding to further compress data size. In this
application example, we mainly focus on sampling to demonstrate the use of MLG
spectrum space to represent structured data.

We test MLG sampling on the RGB icon dataset [17]. This dataset contains
several icon images of size $\mathbb{R}^{16\times 16\times 3}$. Our goal is to
sample a subset of $K=16\times 16\times 3=768$ transformed data points in
spectral domain. We follow Algorithm 1 to sample the transformed data before
using the sampled data to recover the original images. In this case, a 3-layer
MLG is constructed where the intralayer connections are formed by connecting
to the 4-neighborhood pixels and the interlayer connections connect the
counterparts of the same pixel in all layers. We compare our methods with two
GSP-based methods. The first method (GFT) builds a graph with $16\times 16$
nodes and implements GFT on signals in each layer separately. The second
method (GFT2) builds a graph of size $16\times 16$ for pixels that are
neighbor-4 connected and constructs a graph of size $3$ for (fully connected)
frames. Let $\mathbf{E}$ be the spectrum of frame-wise graph, and $\mathbf{F}$
be the spectrum of pixel-wise graph. We apply GFT on the 2-D signal
$\mathbf{s}\in\mathbb{R}^{3\times 256}$ to generate
$\mathbf{\hat{s}}=\mathbf{E}^{\mathrm{T}}\mathbf{s}\mathbf{F}.$ For both GSP
and M-GSP methods, we apply the Laplacian matrix/tensor. We order the
transformed signal in descending order of the transformed values from top left
to bottom right, i.e., keeping larger values in Fig. 2. As mentioned in
Section III-A2, there are three directions for sampling node values. In layer-
wise sampling, we remove samples layer-by-layer from bottom to up. More
specifically, we remove samples from right to left when sampling each layer.
Similarly, in entity-wise sampling, we remove samples from right to left. When
sampling each column, we move from bottom to top. In block-wise sampling, we
preserve samples in top-left blocks. Since there are only three layers, we
present all results of block-wise sampling layer-by-layer in Fig. 2(c)-2(e).

Defining sampling fraction as the ratio of saved samples to total original
signal samples, we measure errors between the recovered and original images.
The results from are shown in Fig. 2, where MLN-EIG and MLN-HOSVD are from MLG
Fourier and singular space respectively. These results confirm that the
proposed MLG methods provide better performance than GSP.

###  IV-B Edge Detection

Convolutions are widely used for edge detection of images. In traditional
processing, one may first transform RGB to gray-scale images for edge
detection using operators such as Prewitt, Canny, and Sobel [18, 19].
Similarly, geometric convolution methods could also detect edges in gray-scale
images [20]. Unlike traditional (graph) signal processing, M-GSP can model RGB
images by using 3-D windows/blocks and define 3-D convolution kernels as shown
in Fig. 4(b). To construct an MLG for edge detection, we also represent each
RGB color by a layer and connect each pixel in one layer with its counterparts
in other layers for interlayer connections. Each pixel is connected with its
four adjacency neighbors for intralayer connections. The benefit of
constructing such graph is to bypass the computation needed to regenerate
spectrum for each window (block). For a $3\times 3\times 3$ kernel, we
construct a 3-layer MLG with 9 nodes in each layer.

(a) Layer-Wise.

(b) Entity-Wise.

(c) Block-Wise.

Figure 3: Options of Sampling Directions for Tensor Signals.

(a) 2-D.

(b) 3-D.

Figure 4: Example of Convolution Kernels. Figure 5: Results of Edge Detection.

In general, one can first extract smooth image features, before applying edge
detection on the difference between original and smoothed images. To define a
graph localization kernel for smoothing, we start from the single-layer graph.
Let $K$ be the number of nodes and the graph convolution kernel be
$\mathbf{c}\in\mathbb{R}^{K}$ where $c_{i}=1$ for $i=k_{1},\cdots,k_{z}$ with
$z\leq K$; otherwise, $c_{i}=0$. The convolution output for a signal
$\mathbf{s}\in\mathbb{R}^{K}$ is
$\mathbf{s^{\prime}}=\mathbf{V}[(\mathbf{V}^{\mathrm{T}}\mathbf{s})*(\mathbf{V}^{\mathrm{T}}\mathbf{c})]=\sum_{j}(\mathbf{f}_{j}\mathbf{f}_{j}^{\mathrm{T}}\mathbf{s}\sum_{i}[f_{j}]_{k_{i}})$,
where $\mathbf{V}=[\mathbf{f}_{1},\cdots,\mathbf{f}_{K}]$ is the graph
spectrum. Different from a pure shift, the convolution kernel enhances
information around nodes $\\{v_{k_{i}}\\}$. With locally-enhanced gray-scale
images, we can replace centroid of a window by the mean of convolution outputs
inside that window to smooth signals further. By shifting the window over the
entire image, one can obtain a locally-enhanced blurred gray-scale image.

Similarly, since our goal is to smooth signals within the 3-D window in MLG,
we can simply localize the information around the center of the block. Here,
we design two 3-D kernels with size $3\times k\times k$, i.e,
$\mathbf{c}_{1}\in\mathbb{R}^{3\times k^{2}}$ with entry in node $v_{2}$ as 1
and all other entries as zero, and $\mathbf{c}_{2}\in\mathbb{R}^{3\times
k^{2}}$ with entries in node $v_{1},v_{2},v_{3}$ as 1 and all other entries as
zero, indexed as 4(b). We then obtain a smoothed gray-scale image from
convolution output and determine its difference from the original gray-scale
image generated from RGB. A threshold is designed here to detect edges based
on the explicit difference between the smoothed image and its original. MLG
spectra are applied here to implement the MLG convolution. Results from the
MLG edge detection can be visually examined in Fig. 5, compared with results
from several classic detection methods (e.g., Canny, Prewitt, Sobel) and a
GSP-based method. The comparison shows that MLG-based methods yield more
smooth edges and show clearer details.

Figure 6: Results of HSI Segmentation.

###  IV-C Hyperspectral Image Segmentation

Spectral clustering [21] is a widely-used method for unsupervised image
segmentation. In this part, we introduce unsupervised hyperspectral image
(HSI) segmentation based on MLG singular space. To construct a MLG, we first
cluster the multi-spectral frames into $M=10$ clusters as layers and construct
$N=100$ superpixels as entities. The new signals in each layer is the mean of
all frames in that cluster. By constructing intralayer and interlayer
connections using Gaussian distance, we generate an $M$-layer MLG with $N$
nodes in each layer. Implementing HOSVD on the adjacency tensor, we can obtain
the entity-wise (superpixel) singular vectors in descending order of the
entity-wise singular values, i.e.,
$\mathbf{E}_{e}=[\mathbf{e}_{1}\cdots\mathbf{e}_{N}]\in\mathbb{R}^{N\times
N}$. Selecting the first $P$ singular tensors as
$\mathbf{P}_{k}=[\mathbf{e}_{1},\cdots,\mathbf{e}_{P}]\in\mathbb{R}^{N\times
P}$ based on the largest singular gap, we can cluster the rows of
$\mathbf{P}_{k}$ into $Q$ groups based on $k$-means clustering. The superpixel
$i$ is clustered into group $j$ if the $i$-th row of $\mathbf{P}_{k}$ is in
group $j$. Labeling the pixels as the same cluster of its superpixel, we
segment the HSI. Let us compare the M-GSP spectral clustering with $k$-means
and GSP spectra clustering in datasets named Indian Pine, Pavia University and
Salinas. The visualization results of the segmented Salinas HSI are shown as
Fig. 6, where M-GSP illustrates more robust results. In addition to visual
results, we also consider quantitative metrics. We compute the edges of
segmentation results from different methods, and compare them to the ground
truth in Table I. M-GSP continues to show superior performance in all HSIs.

TABLE I: Accuracy of Segmentation Boundaries Data |  $k$-means | GSP | M-GSP
---|---|---|---
IndianP | 0.8257 | 0.8298 | 0.8441
Salinas | 0.9208 | 0.9285 | 0.9409
PaviaU | 0.9070 | 0.9088 | 0.9255

##  V Conclusions

In this work, we introduce image modeling based on multilayer graphs. We
present the fundamentals of M-GSP frequency analysis together with the
spectral operations. We further provide several example applications based on
M-GSP operations. Our test results demonstrate the efficacy and strong future
potentials of applying M-GSP in image processing. Other applications of M-GSP
may include dynamic point clouds and video signals.

## References

  * [1] G. Cheung, E. Magli, Y. Tanaka and M. K. Ng, “Graph spectral image processing,” in Proceedings of the IEEE, vol. 106, no. 5, pp. 907-930, May 2018.
  * [2] A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura and P. Vandergheynst, ”Graph signal processing: overview, challenges, and applications,” in Proceedings of the IEEE, vol. 106, no. 5, pp. 808-828, May 2018.
  * [3] G. Fracastoro, D. Thanou and P. Frossard, “Graph transform optimization with application to image compression,” in IEEE Transactions on Image Processing, vol. 29, pp. 419-432, Aug. 2020.
  * [4] M. Onuki, S. Ono, M. Yamagishi and Y. Tanaka, “Graph signal denoising via trilateral filter on graph spectral domain,” in IEEE Transactions on Signal and Information Processing over Networks, vol. 2, no. 2, pp. 137-148, Jun. 2016.
  * [5] M. D. Domenico, A. Solé-Ribalta, E. Cozzo, M. Kivelä, Y. Moreno, M. A. Porter, S. Gómez, and A. Arenas, “Mathematical formulation of multilayer networks,” Physical Review X, vol. 3, no. 4, p. 041022, Dec. 2013.
  * [6] P. Das and A. Ortega, “Graph-based skeleton data compression,” 2020 IEEE 22nd International Workshop on Multimedia Signal Processing (MMSP), Tampere, Finland, Sep. 2020, pp. 1-6.
  * [7] F. Grassi, A. Loukas, N. Perraudin and B. Ricaud, “A time-vertex signal processing framework: scalable processing and meaningful representations for time-series on graphs,” in IEEE Transactions on Signal Processing, vol. 66, no. 3, pp. 817-829, Feb. 2018
  * [8] J. S. Stanley, E. C. Chi and G. Mishne, “Multiway graph signal processing on tensors: integrative analysis of irregular geometries,” in IEEE Signal Processing Magazine, vol. 37, no. 6, pp. 160-173, Oct. 2020.
  * [9] S. Zhang, Q. Deng, and Z. Ding. “Introducing Graph Signal Processing over Multilayer Networks: Theoretical Foundations and Frequency Analysis”, arXiv:2108.13638, 2022. [available online: https://github.com/zsy93/Introducing-Graph-Signal-Processing-over-Multilayer-Networks/blob/main/ms.pdf]
  * [10] A. Afshar, J. C. Ho, B. Dilkina, I. Perros, E. B. Khalil, L. Xiong, and V. Sunderam, “Cp-ortho: an orthogonal tensor factorization framework for spatio-temporal data,” in Proceedings of the 25th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems, Redondo Beach, CA, USA, Jan. 2017, p. 67.
  * [11] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” SIAM Journal on Matrix Analysis and Applications, vol. 21, no. 4, pp. 1253-1278, Jan. 2000.
  * [12] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” ICLR, Toulon, France, Apr. 2017.
  * [13] 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,” in IEEE Signal Processing Magazine, vol. 30, no. 3, pp. 83-98, May 2013.
  * [14] J. Shi, and J. M. Moura, “Graph signal processing: Modulation, convolution, and sampling,” arXiv preprint arXiv:1912.06762, Dec. 2019.
  * [15] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM Review, vol. 51, no. 3, pp. 455-500, Aug. 2009.
  * [16] S. Chen, R. Varma, A. Sandryhaila and J. Kovačević, “Discrete signal processing on graphs: sampling theory,” in IEEE Transactions on Signal Processing, vol. 63, no. 24, pp. 6510-6523, Dec. 2015.
  * [17] S. Zhang, Z. Ding and S. Cui, “Introducing hypergraph signal processing: theoretical foundation and practical applications,” in IEEE Internet of Things Journal, vol. 7, no. 1, pp. 639-660, Jan. 2020.
  * [18] A. K. Cherri, and M. A. Karim, “Optical symbolic substitution: edge detection using Prewitt, Sobel, and Roberts operators,” Applied Optics, vol. 28, no. 21, pp. 4644-4648, Nov. 1989.
  * [19] L. Ding, and A. Goshtasby, “On the Canny edge detector,” Pattern Recognition, vol. 34, no. 3, pp. 721-725, Mar. 2001.
  * [20] S. Zhang, S. Cui and Z. Ding, “Hypergraph-based image processing,” 2020 IEEE International Conference on Image Processing (ICIP), Abu Dhabi, United Arab Emirates, Oct. 2020, pp. 216-220.
  * [21] U. Von Luxburg, “A tutorial on spectral clustering,” Statistics and computing, vol. 17, no. 4, pp. 395-416, Aug. 2007.