# Pyramid Hierarchical Transformer for Hyperspectral Image Classification

Muhammad Ahmad, Muhammad Hassaan Farooq Butt, Manuel Mazzara, Salvatore
Distifano M. Ahmad is with the Department of Computer Science, National
University of Computer and Emerging Sciences, Islamabad, Chiniot-Faisalabad
Campus, Chiniot 35400, Pakistan, and Dipartimento di Matematica e
Informatica—MIFT, University of Messina, Messina 98121, Italy; (e-mail:
mahmad00@gmail.com)M. H. F. Butt is with the Institute of Artificial
Intelligence, School of Mechanical and Electrical Engineering, Shaoxing
University, Shaoxing 312000, China. (e-mail<EMAIL_ADDRESS>Mazzara
is with the Institute of Software Development and Engineering, Innopolis
University, Innopolis, 420500, Russia. (e-mail: m.mazzara@innopolis.ru)S.
Distefano is with Dipartimento di Matematica e Informatica—MIFT, University of
Messina, Messina 98121, Italy. (e-mail<EMAIL_ADDRESS>
###### Abstract

The traditional Transformer model encounters challenges with variable-length
input sequences, particularly in Hyperspectral Image Classification (HSIC),
leading to efficiency and scalability concerns. To overcome this, we propose a
pyramid-based hierarchical transformer (PyFormer). This innovative approach
organizes input data hierarchically into segments, each representing distinct
abstraction levels, thereby enhancing processing efficiency for lengthy
sequences. At each level, a dedicated transformer module is applied,
effectively capturing both local and global context. Spatial and spectral
information flow within the hierarchy facilitates communication and
abstraction propagation. Integration of outputs from different levels
culminates in the final input representation. Experimental results underscore
the superiority of the proposed method over traditional approaches.
Additionally, the incorporation of disjoint samples augments robustness and
reliability, thereby highlighting the potential of our approach in advancing
HSIC.

The source code is available at https://github.com/mahmad00/PyFormer.

###### Index Terms:

Pyramid Network; Hierarchical Transformer; Hyperspectral Image Classification.

##  I Introduction

Hyperspectral Image Classification (HSIC) is crucial in diverse domains,
including remote sensing [1], earth observation [2], urban planning [3],
agriculture [4], forestry [5], target/object detection [6], mineral
exploration [7], environmental monitoring [8, 9], climate change [10], food
processing [11, 12], bakery products [13], bloodstain identification [14, 15],
and meat processing [16, 17]. The abundance of spectral data in Hyperspectral
Images (HSIs) presents challenges and opportunities for classification [18].
While CNNs [19, 20] and its variants [21, 22, 23], and Transformers [24, 25,
26] have shown success in computer vision tasks, there is a growing interest
in exploring Transformer models for advancing HSI analysis.

The vision and spatial-spectral transformers (SSTs) [27, 28, 29, 30, 31, 32,
33, 34] excel in capturing global contextual information via self-attention
mechanisms, facilitating simultaneous consideration of relationships between
all HSI regions [35]. Unlike CNNs, SSTs demonstrate strong scalability to
high-resolution HSIs, effectively handling large datasets without complex
pooling operations. Their adaptability contributes to widespread applicability
in HSIC. Furthermore, SSTs learn stratified representations directly from raw
pixel values, simplifying the model-building process and often leading to
improved performance [36]. The interpretability of attention maps generated by
SSTs aids in understanding model decision-making, highlighting influential
image regions [37].

Despite their success, SSTs have limitations, for instance, training large
SSTs can be computationally demanding [38, 39]. The self-attention mechanism
introduces quadratic complexity with respect to sequence length, potentially
hindering scalability [40]. Unlike CNNs, which inherently possess translation
invariance, SSTs may struggle to capture spatial relationships invariant to
small translations in the input [41]. Moreover, the tokenization process of
dividing input images into fixed-size patches may not efficiently capture
fine-grained details [42, 43]. The quadratic scaling of self-attention poses
challenges, particularly with long sequences. Furthermore, optimal performance
often requires substantial training data, and training on smaller datasets may
lead to overfitting, limiting effectiveness in scenarios with limited labeled
data [25].

Therefore, this work proposed a pyramid hierarchical SST for HSIC. The
hierarchical structure partitions the input into segmentation, each denoting
varying abstraction levels, organized in a pyramid-like manner. Transformer
modules are applied at each level for multi-level processing, ensuring
efficient capture of local and global context. Information flow occurs both
spatially and spectrally within the hierarchy, fostering communication and
abstraction propagation. Integration of Transformer outputs from different
levels yields the final output maps. In short, the following contributions are
made;

  1. 1.

Hierarchical Segmentation: Input sequences are divided into hierarchical
segments, each representing varying levels of abstraction or granularity.

  2. 2.

Pyramid Organization: These segments adopt a pyramid-like structure, wherein
the lowest level retains detailed information while higher levels convey
increasingly abstract representations.

  3. 3.

Multi-level Processing: Transformer modules are independently applied at each
level of the hierarchy, facilitating efficient capture of both local and
global context.

Figure 1: The Pyramid-type Hierarchical Transformer features a structured
pyramid block comprising two levels. Each level incorporates two convolutional
layers with filter sizes of ($32\times S\times S\times B$) and ($64\times
S\times S\times B$), respectively, followed by a residual connection and
activation function. The hierarchical transformer block, receiving the learned
multi-scale information, consists of two layers and four multi-heads. The
acquired information undergoes flattening and is subsequently subjected to the
ReLU activation function and L2 regularization technique. This regularization
aids in mitigating overfitting by reducing weights, rendering the network less
responsive to minor input variations. Finally, the output layer employs
softmax activation for final maps.

##  II Proposed Methodology

An HSI cube, denoted as $X=\\{x_{i},y_{i}\\}\in\mathcal{R}^{(M\times N\times
B)}$, comprises spectral vectors
$x_{i}=\\{x_{i,1},x_{i,2},x_{i,3},\dots,x_{i,L}\\}$, and $y_{i}$ and
corresponding class labels $x_{i}$. The cube is initially divided into
overlapping 3D patches, each centered at spatial coordinates $(\alpha,\beta)$
and spanning $S\times S$ pixels across $B$ bands. The total count of extracted
patches ($m$) from $X$ is $(M-S+1)\times(N-S+1)$, where a patch
$P_{\alpha,\beta}$ covers spatial dimensions within $\alpha\pm\frac{S-1}{2}$
and $\beta\pm\frac{S-1}{2}$. These patches, along with their central pixel
labels, constitute the training $X_{train}$, validation $X_{val}$, and a test
$X_{test}$ sets, ensuring $|X_{train}|\ll|X_{val}|$ and
$|X_{train}|\ll|X_{test}|$ to prevent sample overlaps and biases. The model’s
overall structure is presented in Figure 1.

Let ($S,S,B^{*}$) denote the input shape, $B^{*}$ where represents the reduced
number of bands. The scaling factor, $Scale=2^{Level}$, is computed. This
scale is utilized to derive the input shape for the pyramid structure layers,
given by
$Input~{}shape=(\frac{S}{Scale},\frac{S}{Scale},\frac{B^{*}}{Scale})$. This
input is fed into a convolutional layer to extract spatial-spectral semantic
features from HSI patches. Each patch, with dimensions ($S\times S\times
B^{*}$), undergoes processing using 3D convolutional layers with kernel sizes
($32\times 1\times S\times S$) and ($64\times 1\times S\times S$), along with
ReLU activation. The activation maps for the spatial-spectral position
$(x,y,z)$ at the i-th feature map and j-th layer are denoted as
$f_{i,j}^{(x,y,z)}$ [19];

$f_{(i,j)}(x,y,z)=\sigma\bigg{(}\sum_{m=1}^{M}\sum_{n=1}^{N}\sum_{p=1}^{P}w_{(i,j,m,n,p)}\times\\\
x_{(m+x-1,n+y-1,p+z-1)+b_{i,j}}\bigg{)}$ (1)

Where $f_{i,j}(x,y,z)$ is the activation value, $\sigma$ is the activation
function, $w_{i,j,m,n,p}$ are the weights of the convolutional kernel,
$x_{m+x-1,n+y-1,p+z-1}$ represents the input patch values, and $b_{i,j}$ is
the bias term. Let $f_{i,j}\in\mathbb{R}^{N\times D}$ denote the input tensor
to the transformer, where $N$ represents the number of patches, and $D$
signifies the dimensionality of each patch after the convolutional process.
This encoding is integrated with the input embeddings, augmenting the model
with spatial arrangement details. The foundational architecture of the
transformer centers around the encoder, which consists of multiple layers
incorporating multimodal attention and feedforward convolutional networks. The
attention mechanism assumes a critical role in enabling the model to capture
intricate relationships between distinct patches. Specifically, for a given
input $v_{i,j}$ each layer within the transformer encoder encompasses;

$H_{att}^{l}=Attention(H^{(l-1)})$ (2)

$H_{FF\\_Conv}^{l}=FF\\_Conv(H_{att}^{l})$ (3)

$H^{l}=H^{(l-1)}+H_{FF\\_Conv}^{l}$ (4)

Equations 2, 3, and 4 define the attention, feedforward convolutional neural
network, and residual connections within the transformer layers, respectively.
Later, $H^{l}$ is flattened as $W=reshape(H^{l},(m\times n\times p,1))$, where
$m,n$ and $p$ represent the batch size, height, and width, respectively.
Subsequently, $L_{2}$ regularization $L_{2}Reg=\lambda|W|^{2}_{2}$ is added to
the loss function during training, with ReLU as the activation function and
$\lambda=0.01$ as the regularization parameter. Finally, a softmax function is
employed to generate the classification maps.

##  III Experimental Results and Discussion

To ensure fairness, maintaining identical experimental conditions is crucial,
including consistent geographical locations for model selection and uniform
sample numbers for each training round in cross-validation. Random sample
selection can introduce variability, potentially causing discrepancies among
models executed at different times. Another common issue in recent literature
is the overlap of training and test samples, leading to biased models with
inflated accuracy. To mitigate this, the proposed method ensures that while
training, validation, and test samples are randomly selected, efforts are made
to prevent any overlap between these sets, thereby reducing biases introduced
by overlapping samples. In our experimental setup, the proposed PyFormer was
assessed using a mini-batch size of 128, the Adam optimizer, and specific
learning parameters i.e., learning rate of 0.0001 and a decay rate of 1e-06,
over 50 epochs. We systematically tested various configurations to
comprehensively evaluate the proposed model. This exploration aimed to
thoroughly understand the model’s performance under diverse training scenarios
and spatial resolutions.

TABLE I: Performance on different train ratios across different datasets. Figure 2 illustrates the geographical maps showcasing the best results attained in these experiments. Datasets | Metrics | Data Split Ratio
---|---|---
5% | 10% | 15% | 20% | 25%
PU | OA | 96.28 | 98.33 | 99.35 | 98.02 | 99.80
AA | 93.81 | 97.08 | 98.83 | 97.35 | 99.63
KA | 95.07 | 97.78 | 99.14 | 97.39 | 99.73
F1 | 94.44 | 97.33 | 99.0 | 96.55 | 99.66
SA | OA | 97.53 | 99.08 | 99.33 | 99.75 | 99.86
AA | 98.37 | 94.42 | 99.71 | 99.80 | 99.83
KA | 97.25 | 98.98 | 99.27 | 99.72 | 99.84
F1 | 98.43 | 99.31 | 99.81 | 99.87 | 100.0
HU | OA | 93.11 | 97.36 | 97.13 | 98.19 | 98.18
AA | 92.11 | 95.97 | 96.46 | 97.28 | 98.26
KA | 92.55 | 97.14 | 96.9 | 98.05 | 98.03
F1 | 86.56 | 90.68 | 90.56 | 91.62 | 92.25

(a) Pavia University

(b) Salinas

(c) University of Houston

Figure 2: The PyFormer model achieves kappa accuracies of 99.73%, 99.84%, and 98.01% on the Disjoint Training, Validation, and Test samples for the Pavia University (PU), Salinas (SA), and University of Houston (UH) datasets, respectively. TABLE II: Performance on different patch sizes across different datasets Datasets | Metrics | Different Patch Size
---|---|---
$2\times 2$ | $4\times 4$ | $6\times 6$ | $8\times 8$ | $10\times 10$
PU | OA | 95.88 | 98.72 | 99.45 | 99.22 | 96.28
AA | 94.62 | 98.61 | 99.0 | 99.25 | 95.71
KA | 94.53 | 98.3 | 99.27 | 98.97 | 95.05
F1 | 95 | 98.55 | 99.0 | 99.33 | 94.88
SA | OA | 94.5 | 98.53 | 99.29 | 98.46 | 99.45
AA | 96.26 | 99.44 | 99.71 | 98.81 | 99.51
KA | 93.88 | 98.36 | 99.21 | 98.28 | 99.38
F1 | 95.87 | 99.5 | 99.68 | 98.93 | 99.68
HU | OA | 88.36 | 88.24 | 97.13 | 98.19 | 90.19
AA | 86.86 | 84.34 | 96.46 | 97.28 | 89.03
KA | 87.41 | 87.26 | 96.9 | 98.05 | 89.03
F1 | 81.68 | 79.33 | 90.56 | 91.62 | 83.81
TABLE III: Performance on different Number of Heads across different datasets Datasets | Metrics | Different Number of Heads
---|---|---
2 | 4 | 6 | 8 | 10
PU | OA | 99.33 | 95.93 | 98.38 | 99.11 | 99.36
AA | 98.83 | 95.99 | 96.89 | 98.42 | 98.61
KA | 99.12 | 94.65 | 97.86 | 98.82 | 99.15
F1 | 98.88 | 95.55 | 97.4 | 98.66 | 98.77
SA | OA | 98.67 | 99.5 | 99.33 | 99.33 | 99.61
AA | 98.13 | 99.61 | 98.98 | 99.56 | 99.66
KA | 98.52 | 99.44 | 99.25 | 99.26 | 99.57
F1 | 97.62 | 99.56 | 99.12 | 99.56 | 99.75
HU | OA | 97.49 | 97.67 | 97.82 | 84.47 | 96.77
AA | 96.8 | 96.57 | 97.32 | 85.92 | 96.39
KA | 97.29 | 97.48 | 97.65 | 83.25 | 96.51
F1 | 96.87 | 97.26 | 97.56 | 84.4 | 96.53
TABLE IV: Performance on different Number of Layers across different datasets Datasets | Metrics | Different Number of Layers
---|---|---
2 | 4 | 6 | 8 | 10
PU | OA | 98.79 | 99.37 | 99.44 | 99.11 | 97.94
AA | 98.22 | 98.83 | 99.19 | 99.0 | 95.43
KA | 98.4 | 99.17 | 99.38 | 99.26 | 97.27
F1 | 98.11 | 99 | 99.33 | 99.22 | 96.33
SA | OA | 99.47 | 99.07 | 98.89 | 99.33 | 99.27
AA | 99.68 | 99.53 | 99.87 | 99.24 | 98.63
KA | 99.42 | 98.86 | 99.65 | 98.76 | 98.07
F1 | 99.81 | 99.68 | 99.87 | 97.31 | 98.87
HU | OA | 97.69 | 98.1 | 97.63 | 97.28 | 97.81
AA | 96.55 | 97.47 | 97.03 | 96.50 | 97.13
KA | 97.5 | 97.95 | 97.44 | 97.06 | 97.64
F1 | 97 | 97.75 | 97.4 | 96.6 | 97.26

Initially, we examine four critical factors impacting the model’s performance:
patch sizes and training samples as shown in Tables I and II, and the Number
of heads and layers in the Transformer model as shown in Tables III and IV.
These factors are pivotal for performance optimization. Ensuring an adequately
sized training set covering diverse spectral signatures and representative
samples from each class is crucial. Increasing the number of labeled training
samples offers the model more diverse examples, enhancing learning and
generalization capabilities as shown in Table I. Sufficient samples enhance
model performance, especially in accommodating variations within classes. An
imbalanced distribution of training samples among classes can bias models
towards dominant classes, compromising generalization. Balancing sample
distribution across classes is crucial to mitigate biases and enhance the
model’s generalization ability across all classes. Moreover, patch size
denotes the spatial extent of input patches, crucial for capturing local
spatial information and contextual relationships within HSI data. The choice
of patch size significantly impacts the model’s ability to capture spatial
details and contextual dependencies. Larger patch sizes enhance the model’s
understanding of global spatial relationships, facilitating improved spatial
feature extraction and incorporation of broader spectral information.
Conversely, smaller patch sizes focus on local details, advantageous for
intricate patterns or objects with distinct spatial characteristics and
spectral variations as shown in Table II.

Tables V and VI provide a summary of OA, AA, and kappa, along with the average
for each class. The best performance in each row is highlighted in bold. For
comparison, we have chosen several SOTA models, Transformer [44],
Spectralformer [45], HiT [46], CSiT [47], and WaveFormer (WF) [25]. The
transformer architecture utilizes a standard ViT [48] that has been adapted
for pixel-wise HSI input, consisting of five encoder blocks. CSiT is evaluated
without a Cross-Spectral Attention Fusion module. All results are obtained
using the specified parameter configurations from the original papers to
facilitate direct comparison.

TABLE V: Pavia University: PyFormer is compared against other SOTA models. All
models are evaluated with training/validation/test samples distributed as
5%/5%/90%, respectively.

Class | SF [45] | ViT [44] | WF [25] | CSiT[47] | HiT [46] | PyFormer
---|---|---|---|---|---|---
Asphalt | 92.67% | 95.67% | 96.21% | 93.84% | 95.21% | 95.82%
Meadows | 92.79% | 88.37% | 99.33% | 95.23% | 92.54% | 99.47%
Gravel | 90.60% | 73.71% | 81.31% | 88.79% | 81.18% | 85.61%
Trees | 98.15% | 98.03% | 96.77% | 96.19% | 97.21% | 98.79%
Painted | 98.28% | 99.01% | 100% | 99.18% | 100% | 100%
Soil | 93.29% | 89.26% | 91.55% | 91.99% | 91.93% | 97.18%
Bitumen | 83.01% | 79.40% | 89.55% | 92.06% | 92.75% | 92.88%
Bricks | 84.50% | 85.54% | 89.34% | 82.25% | 87.59% | 89.97%
Shadows | 99.77% | 99.65% | 96.47% | 99.19% | 99.47% | 95.61%
OA | 92.30% | 89.32% | 95.66% | 93.35% | 91.35% | 96.28%
AA | 88.86% | 87.39% | 93.39% | 90.48% | 85.07% | 93.81%
$\kappa$ | 89.66% | 86.60% | 94.22% | 91.13% | 88.94% | 95.07%

TABLE VI: University of Houston: PyFormer is compared against other SOTA
models. All models are evaluated with training/validation/test samples
distributed as 10%/10%/90%, respectively.

Class | SF [45] | ViT [44] | WF [25] | CSiT [47] | HiT [46] | PyFormer
---|---|---|---|---|---|---
Healthy grass | 93.31% | 90.28% | 98.90% | 93.39% | 97.26% | 98.98%
Stressed grass | 97.81% | 98.21% | 97.60% | 99.54% | 97.29% | 99.63%
Synthetic grass | 100% | 96.90% | 99.82% | 100% | 98.74% | 99.71%
Trees | 100% | 100% | 99.19% | 98.09% | 95.78% | 99.48%
Soil | 98.16% | 96.89% | 99.79% | 97.70% | 98.41% | 100%
Water | 100% | 98.21% | 98.07% | 100% | 91.36% | 96.74%
Residential | 87.83% | 82.82% | 90.64% | 90.96% | 94.60% | 97.35%
Commercial | 85.91% | 79.83% | 97.38% | 89.18% | 91.82% | 96.83%
Road | 75.33% | 76.88% | 97.40% | 90.62% | 92.39% | 97.36%
Highway | 82.52% | 80.31% | 97.75% | 93.22% | 90.61% | 99.18%
Railway | 79.19% | 83.17% | 96.76% | 87.91% | 89.09% | 98.84%
Parking Lot 1 | 72.76% | 69.42% | 97.56% | 83.15% | 94.18% | 95.72
Parking Lot 2 | 79.49% | 63.08% | 65.33% | 84.11% | 82.51% | 78.39%
Tennis Court | 93.90% | 90.36% | 98.83% | 97.21% | 91.55% | 100%
Running Track | 97.50% | 94.40% | 99.05% | 100% | 96.72% | 100%
OA | 88.45% | 86.45% | 96.54% | 93.09% | 93.06% | 97.36%
AA | 87.81% | 86.60% | 95.60% | 92.06% | 86.61% | 95.97%
$\kappa$ | 87.50% | 85.35% | 96.26% | 92.53% | 92.50% | 97.14%

The detailed results of the aforementioned models can be found in Tables V and
VI. In summary, the proposed PyFormer model exhibits outstanding performance,
surpassing SOTA ViT-based models across various evaluation metrics, including
OA, AA, and $\kappa$ coefficient. A comprehensive analysis of the quantitative
results indicates that PyFormer consistently achieves superior performance
across different categories, demonstrating significant improvements in
accuracy, as illustrated in the Tables. Notably, while the performance gaps
are relatively small in the PU dataset due to the abundance of samples, the UH
dataset presents a considerable challenge for modeling. For instance, when
evaluating the challenging UH dataset, PyFormer outperforms the baseline ViT
by more than 7% and exceeds SF by approximately 4%. Moreover, the AA achieved
by PyFormer surpasses that of both ViT and SpectralFormer by margins of around
5%, highlighting the potential effectiveness of spatial-spectral feature
extraction. In comparison with the most recent spatial-spectral Transformer
and CSiT models, PyFormer consistently delivers promising results,
demonstrating its proficiency in both spectral and spectral-spatial feature
extraction tasks. It is noteworthy that while HiT excels in identifying land-
cover classes with spectral-spatial information, PyFormer approaches similar
levels of performance. In conclusion, these findings underscore the robustness
and effectiveness of the PyFormer, particularly in scenarios where the
extraction of spatial-spectral information is crucial, especially considering
the limited availability of training samples.

##  IV Conclusions

This paper introduces PyFormer, a novel approach that leverages the strengths
of Pyramid and Vision Transformer for Hyperspectral Image Classification
(HSIC). By extracting multi-scale spatial-spectral features using Pyramid and
integrating them into a transformer encoder, PyFormer can effectively capture
both local texture patterns and global contextual relationships within a
single, end-to-end trainable model. A key innovation is the incorporation of
Pyramid convolutions within the transformer’s attention mechanism,
facilitating enhanced integration of spectral and structural information.
Extensive experiments demonstrate that PyFormer achieves SOTA performance,
particularly on challenging datasets with limited training data. In addition
to superior classification accuracy, PyFormer exhibits robustness and
generalizability, showing promise for addressing real-world problems. Future
research could explore techniques such as self-supervised pre-training and
network optimizations to further enhance PyFormer’s performance, especially in
scenarios with limited data availability.

## References

  * [1] M. Ahmad, S. Shabbir, S. K. Roy, D. Hong, X. Wu, J. Yao, A. M. Khan, M. Mazzara, S. Distefano, and J. Chanussot, “Hyperspectral image classification—traditional to deep models: A survey for future prospects,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , 2021.
  * [2] V. Lodhi, D. Chakravarty, and P. Mitra, “Hyperspectral imaging for earth observation: Platforms and instruments,” _Journal of the Indian Institute of Science_ , vol. 98, pp. 429–443, 2018.
  * [3] Y. Li, D. Hong, C. Li, J. Yao, and J. Chanussote, “Hd-net: High-resolution decoupled network for building footprint extraction via deeply supervised body and boundary decomposition,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 209, pp. 51–65, 2024.
  * [4] B. Lu, P. D. Dao, J. Liu, Y. He, and J. Shang, “Recent advances of hyperspectral imaging technology and applications in agriculture,” _Remote Sensing_ , vol. 12, no. 16, p. 2659, 2020.
  * [5] T. Adão, J. Hruška, L. Pádua, J. Bessa, E. Peres, R. Morais, and J. J. Sousa, “Hyperspectral imaging: A review on uav-based sensors, data processing and applications for agriculture and forestry,” _Remote sensing_ , vol. 9, no. 11, p. 1110, 2017.
  * [6] C. Li, B. Zhang, D. Hong, J. Yao, and J. Chanussot, “Lrr-net: An interpretable deep unfolding network for hyperspectral anomaly detection,” _IEEE Transactions on Geoscience and Remote Sensing_ , 2023.
  * [7] E. Bedini, “The use of hyperspectral remote sensing for mineral exploration: A review,” _Journal of Hyperspectral Remote Sensing_ , vol. 7, no. 4, pp. 189–211, 2017.
  * [8] C. Weber, R. Aguejdad, X. Briottet, J. Avala, S. Fabre, J. Demuynck, E. Zenou, Y. Deville, M. S. Karoui, F. Z. Benhalouche _et al._ , “Hyperspectral imagery for environmental urban planning,” in _IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, 2018, pp. 1628–1631.
  * [9] M. B. Stuart, A. J. McGonigle, and J. R. Willmott, “Hyperspectral imaging in environmental monitoring: A review of recent developments and technological advances in compact field deployable systems,” _Sensors_ , vol. 19, no. 14, p. 3071, 2019.
  * [10] C. B. Pande and K. N. Moharir, “Application of hyperspectral remote sensing role in precision farming and sustainable agriculture under climate change: A review,” _Climate Change Impacts on Natural Resources, Ecosystems and Agricultural Systems_ , pp. 503–520, 2023.
  * [11] M. H. Khan, Z. Saleem, M. Ahmad, A. Sohaib, H. Ayaz, M. Mazzara, and R. A. Raza, “Hyperspectral imaging-based unsupervised adulterated red chili content transformation for classification: Identification of red chili adulterants,” _Neural Computing and Applications_ , vol. 33, no. 21, pp. 14 507–14 521, 2021.
  * [12] M. H. Khan, Z. Saleem, M. Ahmad, A. Sohaib, H. Ayaz, and M. Mazzara, “Hyperspectral imaging for color adulteration detection in red chili,” _Applied Sciences_ , vol. 10, no. 17, p. 5955, 2020.
  * [13] Z. Saleem, M. H. Khan, M. Ahmad, A. Sohaib, H. Ayaz, and M. Mazzara, “Prediction of microbial spoilage and shelf-life of bakery products through hyperspectral imaging,” _IEEE Access_ , vol. 8, pp. 176 986–176 996, 2020\.
  * [14] M. H. F. Butt, H. Ayaz, M. Ahmad, J. P. Li, and R. Kuleev, “A fast and compact hybrid cnn for hyperspectral imaging-based bloodstain classification,” in _2022 IEEE Congress on Evolutionary Computation (CEC)_. IEEE, 2022, pp. 1–8.
  * [15] M. Zulfiqar, M. Ahmad, A. Sohaib, M. Mazzara, and S. Distefano, “Hyperspectral imaging for bloodstain identification,” _Sensors_ , vol. 21, no. 9, p. 3045, 2021.
  * [16] H. Ayaz, M. Ahmad, M. Mazzara, and A. Sohaib, “Hyperspectral imaging for minced meat classification using nonlinear deep features,” _Applied Sciences_ , vol. 10, no. 21, p. 7783, 2020.
  * [17] H. Ayaz, M. Ahmad, A. Sohaib, M. N. Yasir, M. A. Zaidan, M. Ali, M. H. Khan, and Z. Saleem, “Myoglobin-based classification of minced meat using hyperspectral imaging,” _Applied Sciences_ , vol. 10, no. 19, p. 6862, 2020\.
  * [18] D. Hong, B. Zhang, X. Li, Y. Li, C. Li, J. Yao, N. Yokoya, H. Li, P. Ghamisi, X. Jia, A. Plaza, P. Gamba, J. A. Benediktsson, and J. Chanussot, “Spectralgpt: Spectral remote sensing foundation model,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2024, dOI:10.1109/TPAMI.2024.3362475.
  * [19] M. Ahmad, A. M. Khan, M. Mazzara, S. Distefano, M. Ali, and M. S. Sarfraz, “A fast and compact 3-d cnn for hyperspectral image classification,” _IEEE Geoscience and Remote Sensing Letters_ , 2020.
  * [20] U. Ghous, M. S. Sarfraz, M. Ahmad, C. Li, and D. Hong, “(2+1)d extreme xception net for hyperspectral image classification,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , pp. 1–14, 2024\.
  * [21] M. Ahmad and M. Mazzara, “Scsnet: Sharpened cosine similarity-based neural network for hyperspectral image classification,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 21, pp. 1–4, 2024.
  * [22] D. Hong, L. Gao, J. Yao, B. Zhang, A. Plaza, and J. Chanussot, “Graph convolutional networks for hyperspectral image classification,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 59, no. 7, pp. 5966–5978, 2021.
  * [23] A. Jamali, S. K. Roy, D. Hong, P. M. Atkinson, and P. Ghamisi, “Attention graph convolutional network for disjoint hyperspectral image classification,” _IEEE Geoscience and Remote Sensing Letters_ , pp. 1–1, 2024.
  * [24] J. Yao, B. Zhang, C. Li, D. Hong, and J. Chanussot, “Extended vision transformer (exvit) for land use and land cover classification: A multimodal deep learning framework,” _IEEE Transactions on Geoscience and Remote Sensing_ , 2023.
  * [25] M. Ahmad, U. Ghous, M. Usama, and M. Mazzara, “Waveformer: Spectral–spatial wavelet transformer for hyperspectral image classification,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 21, pp. 1–5, 2024.
  * [26] X. Huang, M. Dong, J. Li, and X. Guo, “A 3-d-swin transformer-based hierarchical contrastive learning method for hyperspectral image classification,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 60, pp. 1–15, 2022.
  * [27] L. Wang, Z. Zheng, N. Kumar, C. Wang, F. Guo, and P. Zhang, “Multilevel class token transformer with cross tokenmixer for hyperspectral images classification,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 62, pp. 1–13, 2024.
  * [28] X. He, Y. Chen, and Z. Lin, “Spatial-spectral transformer for hyperspectral image classification,” _Remote Sensing_ , vol. 13, no. 3, 2021.
  * [29] G. Wang, Y. Wang, Z. Pan, X. Wang, J. Zhang, and J. Pan, “Vitfsl-baseline: A simple baseline of vision transformer network for few-shot image classification,” _IEEE Access_ , vol. 12, pp. 11 836–11 849, 2024.
  * [30] Y. Ma, Y. Lan, Y. Xie, L. Yu, C. Chen, Y. Wu, and X. Dai, “A spatial-spectral transformer for hyperspectral image classification based on global dependencies of multi-scale features,” _Remote Sensing_ , vol. 16, no. 2, 2024.
  * [31] J. Lian, L. Wang, H. Sun, and H. Huang, “Gt-had: Gated transformer for hyperspectral anomaly detection,” _IEEE Transactions on Neural Networks and Learning Systems_ , pp. 1–15, 2024.
  * [32] S. Mei, Z. Han, M. Ma, F. Xu, and X. Li, “A novel center-boundary metric loss to learn discriminative features for hyperspectral image classification,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 62, pp. 1–16, 2024.
  * [33] A. Jamali, S. K. Roy, D. Hong, P. M. Atkinson, and P. Ghamisi, “Spatial-gated multilayer perceptron for land use and land cover mapping,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 21, pp. 1–5, 2024.
  * [34] Y. Xiao, Q. Yuan, K. Jiang, J. He, C.-W. Lin, and L. Zhang, “Ttst: A top-k token selective transformer for remote sensing image super-resolution,” _IEEE Transactions on Image Processing_ , vol. 33, pp. 738–752, 2024.
  * [35] S. Mei, C. Song, M. Ma, and F. Xu, “Hyperspectral image classification using group-aware hierarchical transformer,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 60, pp. 1–14, 2022.
  * [36] D. Hong, Z. Han, J. Yao, L. Gao, B. Zhang, A. Plaza, and J. Chanussot, “Spectralformer: Rethinking hyperspectral image classification with transformers,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 60, pp. 1–15, 2022.
  * [37] J. Chen, C. Yang, L. Zhang, L. Yang, L. Bian, Z. Luo, and J. Wang, “Tccu-net: Transformer and cnn collaborative unmixing network for hyperspectral image,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , pp. 1–20, 2024.
  * [38] J. Fang, J. Yang, A. Khader, and L. Xiao, “Mimo-sst: Multi-input multi-output spatial-spectral transformer for hyperspectral and multispectral image fusion,” _IEEE Transactions on Geoscience and Remote Sensing_ , pp. 1–1, 2024.
  * [39] M. Ye, J. Chen, F. Xiong, and Y. Qian, “Adaptive graph modeling with self-training for heterogeneous cross-scene hyperspectral image classification,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 62, pp. 1–15, 2024.
  * [40] W. Huang, Y. Deng, S. Hui, Y. Wu, S. Zhou, and J. Wang, “Sparse self-attention transformer for image inpainting,” _Pattern Recognition_ , vol. 145, p. 109897, 2024.
  * [41] Y. Sun, X. Zhi, S. Jiang, G. Fan, X. Yan, and W. Zhang, “Image fusion for the novelty rotating synthetic aperture system based on vision transformer,” _Information Fusion_ , vol. 104, p. 102163, 2024.
  * [42] T. Kim, J. Kim, H. Oh, and J. Kang, “Deep transformer based video inpainting using fast fourier tokenization,” _IEEE Access_ , vol. 12, pp. 21 723–21 736, 2024.
  * [43] Y. Shi, J. Xia, M. Zhou, and Z. Cao, “A dual-feature-based adaptive shared transformer network for image captioning,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 73, pp. 1–13, 2024.
  * [44] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [45] D. Hong, Z. Han, J. Yao, L. Gao, B. Zhang, A. Plaza, and J. Chanussot, “Spectralformer: Rethinking hyperspectral image classification with transformers,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 60, pp. 1–15, 2021.
  * [46] X. Yang, W. Cao, Y. Lu, and Y. Zhou, “Hyperspectral image transformer classification networks,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 60, pp. 1–15, 2022.
  * [47] W. He, W. Huang, S. Liao, Z. Xu, and J. Yan, “Csit: A multiscale vision transformer for hyperspectral image classification,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 15, pp. 9266–9277, 2022.
  * [48] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, J. Uszkoreit, and N. Houlsby, “An image is worth 16x16 words: Transformers for image recognition at scale,” 2021.


11institutetext: University of Trieste, Italy 22institutetext: University of
Pisa, Italy

# Retrieval-Augmented Mining of Temporal Logic Specifications from Data

Gaia Saveri 11 2 2 Luca Bortolussi 11

###### Abstract

The integration of cyber-physical systems (CPS) into everyday life raises the
critical necessity of ensuring their safety and reliability. An important step
in this direction is requirement mining, i.e. inferring formally specified
system properties from observed behaviors, in order to discover knowledge
about the system. Signal Temporal Logic (STL) offers a concise yet expressive
language for specifying requirements, particularly suited for CPS, where
behaviors are typically represented as time series data. This work addresses
the task of learning STL requirements from observed behaviors in a data-driven
manner, focusing on binary classification, i.e. on inferring properties of the
system which are able to discriminate between regular and anomalous behaviour,
and that can be used both as classifiers and as monitors of the compliance of
the CPS to desirable specifications. We present a novel framework that
combines Bayesian Optimization (BO) and Information Retrieval (IR) techniques
to simultaneously learn both the structure and the parameters of STL formulae,
without restrictions on the STL grammar. Specifically, we propose a framework
that leverages a dense vector database containing semantic-preserving
continuous representations of millions of formulae, queried for facilitating
the mining of requirements inside a BO loop. We demonstrate the effectiveness
of our approach in several signal classification applications, showing its
ability to extract interpretable insights from system executions and advance
the state-of-the-art in requirement mining for CPS.

###### Keywords:

Time-series data Requirement Mining Temporal Logic

##  1 Introduction

The adoption of cyber-physical systems (CPS) is pervasive nowadays: wearable
devices, internet of things applications and domotics systems are part of our
daily lives. However, formal specifications describing their expected
behaviour are not always (completely) available. This is extremely serious
when CPS are deployed in safety-critical scenarios, such as autonomous
vehicles or medical devices, where safety and reliability are key desiderata.
Requirement (or specification) mining is the the task of inferring system
properties from observations of its behaviour [7]. In the context of CPS, such
behaviours are often recorded as time series describing the dynamics of the
system over time (e.g. readings from a sensor at regular time intervals),
hence Signal Temporal Logic (STL) [16, 23] is a popular language for
expressing requirements [6]. STL is indeed a formalism which allows to reason
over real-valued trajectories in a concise yet rich way: for example in STL
one can state properties like "the temperature of the room will reach $25$
degrees within the next $10$ minutes and will stay above $22$ degrees for the
next hour". In this work, we tackle the objective of learning STL requirements
from observed behaviours of a system, in a data-driven approach: differently
from formal methods algorithms, our proposed solution does not require a
complete model of the system, but only a fixed dataset of recorded
trajectories. Specifically, we frame the STL requirement mining problem in the
binary classification scenario, i.e. our goal is to learn a STL formula which
is able to discriminate between regular and anomalous trajectories. In this
context the mined STL specification can be used both as a classifier and as a
monitor for the system at hand, and leveraged for maintenance and bug
detection as well as a tool for system modeling and comprehension.

##### Our contribution

consists in proposing a novel framework for STL specification mining, which
involves combining Bayesian Optimization (BO) and Information Retrieval (IR)
techniques, to simultaneously learn the structure and the parameters of the
specifications, without neither the need of a (generative) model of the system
under analysis, nor of restrictions on the STL grammar to some specific
fragment. Indeed we: (i) leverage the semantics of STL for mapping formulae in
a continuous semantic-preserving latent space, using the method of [9], see
Section 2; (ii) build a dense vector database containing the embeddings of
millions of formulae, carefully constructed in an incremental way, as detailed
in Section 3.1. Notably, this step is done once, and the same database is
queried in all the application cases; (iii) frame the requirement mining task
as an optimization problem in the semantic space of STL formulae, using
queries to the database constructed at the previous step for inverting the
embeddings, i.e. mapping continuous semantic vectors back to formulae
(Sections 3). We show the effectiveness of the proposed approach in several
binary classification applications, proving that it is able to extract
interpretable information from system executions, offering insights on the
system under analysis.

##  2 Background

##### Signal Temporal Logic

(STL) is a linear-time temporal logic which expresses properties on
trajectories over dense time intervals [16]. We define as trajectories the
functions $\xi:I\rightarrow D$, where $I\subseteq\mathbb{R}_{\geq 0}$ is the
time domain and $D\subseteq\mathbb{R}^{k},k\in\mathbb{N}$ is the state space.
The syntax of STL is given by:

$\varphi:=tt\mid\pi\mid\lnot\varphi\mid\varphi_{1}\land\varphi_{2}\mid\varphi_{1}\mathbf{U}_{[a,b]}\varphi_{2}$

where $tt$ is the Boolean _true_ constant; $\pi$ is an _atomic predicate_ ,
i.e. a function over variables $\bm{x}\in\mathbb{R}^{n}$ of the form
$f_{\pi}(\bm{x})\geq$0$$ (we refer to $n$ as the number of variables of a STL
formula); $\lnot$ and $\land$ are the Boolean _negation_ and _conjunction_ ,
respectively (from which the _disjunction_ $\lor$ follows by De Morgan’s law);
$\mathbf{U}_{[a,b]}$, with $a,b\in\mathbb{Q},a<b$, is the _until_ operator,
from which the _eventually_ $\mathbf{F}_{[a,b]}$ and the _always_
$\mathbf{G}_{[a,b]}$ temporal operators can be deduced. We can intuitively
interpret the temporal operators over $[a,b]$ as follows: a property is
eventually satisfied if it is satisfied at some point inside the temporal
interval, while a property is globally satisfied if it is true continuously in
$[a,b]$; finally the until operator captures the relationship between two
conditions $\varphi,\psi$ in which the first condition $\varphi$ holds until,
at some point in $[a,b]$, the second condition $\psi$ becomes true. We call
$\mathcal{P}$ the set of well-formed STL formulae. STL is endowed with both a
_qualitative_ (or Boolean) semantics, giving the classical notion of
satisfaction of a property over a trajectory, i.e. $s(\varphi,\xi,t)=$1$$ if
the trajectory $\xi$ at time $t$ satisfies the STL formula $\varphi$, and a
_quantitative_ semantics, denoted by $\rho(\varphi,\xi,t)$. The latter, also
called _robustness_ , is a measure of how robust is the satisfaction of
$\varphi$ w.r.t. perturbations of the signals. Robustness is recursively
defined as:

$\displaystyle\rho(\pi,\xi,t)$ $\displaystyle=$ $\displaystyle
f_{\pi}(\xi(t))\qquad\text{for }\pi(\bm{x})=\big{(}f_{\pi}(\bm{x})\geq
0\big{)}$ $\displaystyle\rho(\lnot\varphi,\xi,t)$ $\displaystyle=$
$\displaystyle-\rho(\varphi,\xi,t)$
$\displaystyle\rho(\varphi_{1}\land\varphi_{2},\xi,t)$ $\displaystyle=$
$\displaystyle\min\big{(}\rho(\varphi_{1},\xi,t),\rho(\varphi_{2},\xi,t)\big{)}$
$\displaystyle\rho(\varphi_{1}\mathbf{U}_{[a,b]}\varphi_{2},\xi,t)$
$\displaystyle=$
$\displaystyle\max_{{t^{\prime}\in[t+a,t+b]}}\big{(}\min\big{(}\rho(\varphi_{2},\xi,t^{\prime}),\min_{{t^{\prime\prime}\in[t,t^{\prime}]}}\rho(\varphi_{1},\xi,t^{\prime\prime})\big{)}\big{)}$

Since they will be used later, we report also the definition of robustness of
derived temporal operators: eventually
$\rho(\mathbf{F}_{[a,b]}\varphi,\xi,t)=\max_{t^{\prime}\in[t+a,t+b]}\rho(\varphi,\xi,t)$
and globally
$\rho(\mathbf{G}_{[a,b]}\varphi,\xi,t)=\min_{t^{\prime}\in[t+a,t+b]}\rho(\varphi,\xi,t)$.
Robustness is compatible with satisfaction via the following _soundness_
property: if $\rho(\varphi,\xi,t)>0$ then $s(\varphi,\xi,t)=$1$$ and if
$\rho(\varphi,\xi,t)<$0$$ then $s(\varphi,\xi,t)=$0$$. When
$\rho(\varphi,\xi,t)=$0$$ arbitrary small perturbations of the signal might
lead to changes in satisfaction value. We omit $t$ from the previous notations
when properties are evaluated at time $t=$0$$.

In this context stochastic processes are probability spaces defined as
triplets $\bm{X}=(\mathcal{T},\mathcal{A},\mu)$ of a trajectory space
$\mathcal{T}$ and a probability measure $\mu$ on a $\sigma$-algebra
$\mathcal{A}$ over $\mathcal{T}$. Given a STL formula $\varphi$ with
predicates interpreted over state variables of $\bm{X}$ and a trajectory of
the stochastic system $\xi(t)$, its robustness $\rho(\varphi,\xi)$ is a
functional $R_{\varphi}:\mathcal{T}\rightarrow\mathbb{R}$, which defines the
real-valued random variable $R_{\varphi}=R_{\varphi}(\bm{X})$, following the
distribution:

$\mathbb{P}(R_{\varphi}(\bm{X})\in[a,b])=\mathbb{P}(\bm{X}\in\\{\xi\in\mathcal{T}|\rho(\varphi,\xi)\in[a,b]\\})$
(1)

The expected value of the previous, i.e. $\mathbb{E}[R_{\varphi}|\bm{X}]$,
intuitively gives a measure of how strongly a formula is satisfied by
trajectories sampled from the process $\bm{X}$ (the higher the value the more
robust the satisfaction) [5]; in the remainder of this paper, this expectation
will be approximated by Monte Carlo sampling.

##### Finding continuous representations

of STL formulae is performed with an ad-hoc kernel in [9], by leveraging the
quantitative semantics of STL. Indeed, robustness allows formulae to be
considered as functionals mapping trajectories into real numbers, i.e.
$\rho(\varphi,\cdot):\mathcal{T}\rightarrow\mathbb{R}$ such that
$\xi\mapsto\rho(\varphi,\xi)$. Considering these as feature maps, and fixing a
probability measure $\mu_{0}$ on the space of trajectories $\mathcal{T}$, a
kernel function capturing similarity among STL formulae on mentioned feature
representations can be defined as:

$k(\varphi,\psi)=\langle\rho(\varphi,\cdot),\rho(\psi,\cdot)\rangle=\int_{\xi\in\mathcal{T}}\rho(\varphi,\xi)\rho(\psi,\xi)d\mu_{0}(\xi)$
(2)

opening the doors to the use of the scalar product in the Hilbert space
$L^{2}$ as a kernel for $\mathcal{P}$; intuitively this results in a kernel
having high positive value for formulae that behave similarly on high-
probability trajectories (w.r.t. $\mu_{0}$), and viceversa low negative value
for formulae that on those trajectories disagree.

Algorithm 1 Sampling a trajectory over the interval $[a,b]$ according to
$\mu_{0}$

$\Delta$, $a$, $b$, $m^{\prime}$, $m^{\prime\prime}$, $\sigma^{\prime}$,
$\sigma^{\prime\prime}$, $q$

$\xi$

$\triangleright$ sample the starting point

$\xi_{0}\sim\mathcal{N}(m^{\prime},\sigma^{\prime})$

$\xi(t_{0})\leftarrow\xi_{0}$

$\triangleright$ sample the total variation

$K\sim(\mathcal{N}(m^{\prime\prime},\sigma^{\prime\prime}))^{2}$

$y_{1},...,y_{N-1}\sim\mathbb{U}([0,K])$

$y_{0}\leftarrow 0$, $y_{n}\leftarrow K$

orderAndRename($y_{0},\ldots,y_{n}$)

$\triangleright$ now $y_{1}\leq y_{2}\leq...\leq y_{N-1}$

$s_{0}\sim\text{Discr}(-1,1)$

while $i\leq N$ do

$s\leftarrow$Binomial($q$) $\triangleright$ $P(s=-1)=q$

$s_{i+1}=s_{i}\cdot s$

$\xi(t_{i+1})=\xi(t_{i})+s_{i+1}(y_{i+1}-y_{i})$

end while

Intuitively, $\mu_{0}$ makes _simple_ trajectories more probable, considering
total variation and number of changes in monotonicity as indicators of
complexity of signals. The measure $\mu_{0}$, operating on piece-wise linear
functions over the interval $\mathcal{I}=[a,b]$ (which is a dense subset of
the set of continuous functions over $\mathcal{I}$), can be algorithmically
defined as in Algorithm 1 (default parameters are set as
$a=0,b=100,\Delta=1,m^{\prime}=m^{\prime\prime}=0^{\prime\prime},\sigma^{\prime}=\sigma^{\prime\prime}=1,q=0.1$).
Note that, although the feature space $\mathbb{R}^{\mathcal{T}}$ (which we
call the _latent semantic space_) into which $\rho$ (and thus Equation (2))
maps formulae is infinite-dimensional, in practice the kernel trick allows to
circumvent this issue by mapping each formula to a vector of dimension equal
to the number of formulae which are in the training set used to evaluate the
kernel (Gram) matrix.

##### Bayesian Optimization (BO)

is a widely-used framework for black-box optimization problems, e.g. of the
form $\bm{x}^{\star}=\text{argmax}_{\bm{x}\in\mathcal{X}}f(\bm{x})$ being
$\mathcal{X}\subset\mathbb{R}^{d}$, which consists in sequentially at each
iteration $t$ observing noisy evaluations of the objective function
$y_{t}=f(\bm{x}_{t})+\varepsilon_{t}$ where the noise is
$\varepsilon_{t}\sim\mathcal{N}(0,\sigma^{2})$, at candidate points
$\bm{x}_{t}$ chosen by an acquisition function. The Gaussian Process Upper
Confidence Bound (GP-UCB) algorithm consists in interpolating the available
observations using a Gaussian Process (GP, i.e. a Bayesian non-parametric
regression approaches) [24], which is the so-called emulation phase, and
adding candidate points to the training set of the GP in an iterative fashion
using Upper-Confidence Bound (UCB, [26]) as acquisition function. More
specifically it considers $f$ to be a sample path from a GP with $0$ mean and
stationary kernel function $k:\mathcal{X}\times\mathcal{X}\mapsto\mathbb{R}$,
i.e. $f\sim GP(0,k)$, so that at each iteration the posterior mean and
variance of $f(\bm{x})$ are derived, and respectively denoted as
$\mu_{t-1}(\bm{x})$ and $\sigma^{2}_{t-1}(\bm{x})$. The UCB acquisition
function selects the next candidate point to be evaluated as
$\bm{x}_{t}=\text{argmax}_{x\in\mathcal{X}}\mu_{t-1}(\bm{x})+\beta_{t}^{\frac{1}{2}}\sigma_{t-1}(\bm{x})$,
with $\beta_{t}\in\mathbb{R}$.

##  3 Retrieval-Augmented STL Requirement Mining

The problem we address with our proposed methodology is that of mining a STL
specification from a dataset of labeled trajectories (where labels indicate
whether the system complies to some desired behaviour), without the need of
inferring or knowing a generative model of the system. This is a supervised
two-class classification problem, in which the input is a set of trajectories
partitioned in the subset of those labeled as regular (or positive) and those
which are anomalous (or negative), denoted respectively as $\mathcal{D}_{p}$
and $\mathcal{D}_{n}$, and the output is a (set of) STL formula(e) able to
discriminate among the two sets. We also assume that such datasets come from
unknown stochastic processes $\bm{X}_{p}$ and $\bm{X}_{n}$, respectively. We
adopt the approach of [19] and optimize the following function in order to
mine a STL formula $\varphi$ with the above desiderata:

$G(\varphi)=\frac{\mathbb{E}[R_{\varphi}|\bm{X}_{p}]-\mathbb{E}[R_{\varphi}|\bm{X}_{n}]}{\sigma(R_{\varphi}|\bm{X}_{p})+\sigma(R_{\varphi}|\bm{X}_{n})}$
(3)

with $R_{\varphi}$ defined in Equation 1. The intuition behind Equation 3 is
that the higher the value of $G(\varphi)$ (i.e. the bigger the difference
between the average robustness of the two stochastic processes, relative to
the length scale given by their standard deviation $\sigma$), the more
robustly the inferred $\varphi$ will separate the sets of positive and
negative trajectories. The key insights of our approach are that: (i) we frame
the learning problem as the optimization of Equation 3 in the latent semantic
space of formulae, i.e. in the space of embeddings of formulae individuated by
the STL kernel of Equation 2; (ii) we construct a dense (hierarchical) vector
database containing the embeddings of millions of STL formulae, that we query
anytime we need to invert the embeddings resulting from the previous step.
This IR step approximates the inverse of Equation 2 by finding specification
whose latent representation is the nearest to the searched embedding (among
the ones in the database), to obtain STL formulae for evaluating $G(\varphi)$.
Hence the focus of our procedure is on inferring the semantic features a
formula needs to have to discriminate between positive and negative
trajectories, without a-priori constraints on its syntactic structure.

Summarizing, the proposed methodology consists in interleaving steps of
Bayesian Optimization (BO) and Information Retrieval (IR): the details of
these two building blocks are given in Sections 3.2 and 3.1, respectively.

###  3.1 Building and Querying the Semantic Vector Database

The STL kernel of Equation 2 allows to map formulae from the discrete space of
their syntactic representation to a continuous space preserving their
semantics (i.e. in the individuated latent space semantically similar formulae
are mapped to nearby representations), however such map is not invertible.
Given the recent advances in indexing structures and search mechanisms for
dense vector databases [12], often accompanied by open-source software, we
devise a simple yet effective strategy for inverting the STL kernel
representations: storing in a reasoned and hierarchical database the
embeddings of millions of formulae (more precisely of pairs containing both
the formula and the corresponding representation), that we call semantic
vector database, and using Approximate Nearest Neighbors (ANN) search
techniques to obtain STL requirements starting from a candidate embedding.

##### Selecting STL Formulae

In order to build such database we first need to select the set of formulae to
store: we do so by fixing the maximum number $M$ of nodes and the maximum
number of variables $N$ (i.e. the maximum dimensionality of the signals)
allowed, then we enumerate all templates (i.e. STL formulae where constants
are replaced by parameters) satisfying those constraints as detailed in
Algorithm 2.

Algorithm 2 Algorithms for generating STL formulae templates $M$, $N$ all_phis
$\triangleright$ STL templates of formulae with max $N$ vars and $M$ nodes
$\text{all\\_phis}\leftarrow[]$ all_phis.append(generateAtomicPropositions())
$\triangleright$ $x_{i}\leq 0$ or $x_{i}\geq 0$, $\forall i\leq N$ for $2\leq
m\leq M$ do $\text{prev\\_phis}\leftarrow$getPhisGivenNodes($m-1$)
$\triangleright$ retrieve templates with $m-1$ nodes
$\text{unary\\_ops}\leftarrow$expandByUnaryOperators(prev_phis)
$\triangleright$ $F$, $G$, $\neg$ all_phis.append(unary_ops)  l_list, r_list
$\leftarrow$getPairsGivenSum($m$) $\triangleright$ all pairs
$(l,r):l+r=m,l\leq r$ for $(l,r)\in$ [l_list, r_list] do
l_phis$\leftarrow$getPhisGivenNodes($l$)
r_phis$\leftarrow$getPhisGivenNodes($r$)
binary_ops$\leftarrow$expandByBinaryOpeators(l_phis, r_phis) $\triangleright$
$\wedge$, $\vee$, $U$ all_phis.append(binary_ops)  end for end for

Once all the templates have been generated, we instantiate each of them on a
set of parameter configurations $\mathcal{P}=\\{p_{0},\ldots
p_{|\mathcal{P}|}\\}$ to be filtered by Signature-based Optimization as
defined in [18]. Denoting as $\varphi(p)$ the application of the parameter set
$p$ to the template formula $\varphi$, and given a set of trajectories
$\hat{\mathcal{T}}=\\{\xi_{0},\ldots,\xi_{|\hat{\mathcal{T}}|}\\}$ (with
$p\in\mathcal{P}$, and $\hat{\mathcal{T}}\subset\mathcal{T}$), it consists in
building a matrix $S\in\mathbb{R}^{|\mathcal{P}|\times|\hat{\mathcal{T}}|}$,
whose rows are called the signature of the formula, s.t.
$S(i,j)=\rho(\varphi(p_{i}),\xi_{j})$, i.e. the cell $(i,j)$ of such matrix
holds the robustness of the $i$-th instantiation of $\varphi$ on the $j$-th
trajectory. In such matrix, the distance between rows can be considered as a
proxy for semantic similarity among different parametrizations of the same
template, hence a selection on formulae can be done by keeping only those
whose signature has a distance higher than a certain threshold
$\tau\in\mathbb{R}_{\geq 0}$ from formulae already selected. This helps in
removing redundancy in the final database structure.

##### Computing Formulae Embeddings

Once the list of formulae to be stored in the database is provided, the STL
kernel embeddings of such specifications needs to be computed, as for Equation
2. In order to do so, one needs to (i) fix a probability measure on the space
of trajectories and (ii) sample a set of $n_{\text{train}}$ formulae against
which the kernel is computed. For the first point, we use the default measure
$\mu_{0}$, computed as in Algorithm 1, proved to be quite general in [9]. For
the second point we instead sample from a distribution $\mathcal{F}$ over STL
formulae which is algorithmically defined by a syntax-tree random recursive
growing scheme, that recursively generates the nodes of a formula given the
probability $p_{\mathit{leaf}}$ of each node being an atomic predicate, and a
uniform distribution over the other operator nodes. We remark that the
$n_{\text{train}}$ determines, by construction, the dimension of the vectors
stored in the database.

##### Index Structure Organization

It is worth noting that the semantic vector database is built once (although
it can be further extended) and leveraged in all the applications, hence it is
fundamental to organize it in a reasoned way, given the different contexts it
can be deployed in. In particular, we hierarchically index the embeddings
based on the number of different variables and the number of nodes appearing
in the formula. The rationale behind the first choice is that not all signals
that we wish to classify with our methodology have the same dimension, the
second choice moves instead from the point of view of interpretability: the
smaller the syntax tree of a formula, the higher its human-understandability,
hence we incorporate in the retrieval strategy the preference for
syntactically simpler formulae.

##### Semantic Database Implementation

For indexing the embeddings, we use the tools provided by the FAISS Python
library [13], which allows us to leverage GPU acceleration for both indexing
and querying steps. Since the amount of vectors we need to store is of the
order of millions, we opt for the Inverted File Index with Product
Quantization (IVFPQ) to accelerate the search procedure (numerical results
will follow in the remainder of this Section), which is a composite index that
both partitions the embeddings in Voronoi cells and reduces their dimension
using Product Quantization [12]. Euclidean (or $L_{2}$) distance is used as
similarity measure between vectors. We remark that the framework we propose is
flexible w.r.t. the choice of the index structure, hence this component can be
replaced with any similarity search engine on vector databases without
altering the overall structure of the algorithm.

Table 1: Quantiles of $AP@K$, $NDCG@K$ and kernel similarity
$k(\varphi,\hat{\varphi})$ between the searched formula $\varphi$ and the top
result $\hat{\varphi}$, across $1000$ queries to the semantic database,
varying the number of nodes in the syntactic tree of the searched formulae.

$AP@K$ | $NDCG@K$ | $k(\varphi,\hat{\varphi})$
---|---|---
$n_{\text{nodes}}$ | 1quart | median | 3quart | 99perc | 1quart | median | 3quart | 99perc | 1quart | median | 3quart | 99perc
$(0,5]$ | $0.804$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $0.810$ | $0.901$ | $0.976$ | $0.992$
$(5,10]$ | $0.763$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $0.774$ | $0.883$ | $0.946$ | $0.979$
$(10,15]$ | $0.684$ | $0.932$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $0.728$ | $0.854$ | $0.928$ | $0.972$
$(15,20]$ | $0.635$ | $0.873$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $1.000$ | $0.742$ | $0.848$ | $0.917$ | $0.976$

In order to build the index we enumerate all formulae templates following
Algorithm 2 with $M=5$ and $N=3$. Given that the measure $\mu_{0}$ is used for
computing the embeddings, we keep into account the typical range of its
trajectories (see Figure 2) for setting the variable thresholds, hence using
the grid of $10$ equally-spaced points between $-4$ and $4$ as possible
parameter instances; for the temporal operators thresholds, we instead
leverage the fact that signals for computing the STL kernel representations
have $100$ points and use the grid of $10$ equally-spaced points between $0$
and $100$ as possible instances for both left and right time-bounds. We select
the formulae to store in the index by signature-based optimization with cosine
distance and threshold $\tau=0.9$. Finally, as done in [9], we use
$n_{\text{train}}=1000$ STL formulae for computing the embeddings. This leads
to a total of $\sim$10\,000\,000$$ stored formulae, that we organize in a
total of $6$ indexes: first splitting them based on the number of variables
they contain (namely $1$, $2$ or $3$), then splitting each of them based on
the maximum number of nodes allowed (namely $4$ or $5$). Of course, alongside
the embeddings, we store the STL formulae organized in a list structure
mirroring the organization of the corresponding representations. Searching for
$1000$ formulae with $3$ variables on such index (hence scanning all the
database) requires $1\pm 0.5s$, while computing their embeddings takes $6\pm
2s$ exploiting GPU acceleration, finally the whole process of creating the
database takes $\sim 20h$ on a AMD EPYC $7542$ machine with $64$ GB of RAM and
a NVidia A$100$ with $80$GB of memory. Moreover, storing the index requires
$\sim 180MB$, which is order of magnitudes less then the $\sim 70GB$ required
in case plain embeddings should be stored (e.g. for exact search).

##### Semantic Database Effectiveness

In order to assess the effectiveness of our information retrieval system, we
use Average Precision at $K$ ($AP@K$) and Normalized Discounted Cumulative
Gain at $K$ ($NDCG@K$), to measure both the relevancy of the results and the
ranking quality. Defining the precision at $K$ (denoted as $P@K$) as the
number of relevant results among the top $K$ retrieved, and $rel_{k}$ as the
indicator function being $1$ when the $k$-th ranked result is relevant and $0$
otherwise, then $AP@K=\frac{\sum_{k=1}^{K}(P@K\cdot
rel_{k})}{\sum_{k=1}^{K}rel_{k}}$. If we instead define $r_{k}$ as the graded
relevance (or gain) of the result at position $k$, i.e. an integer score
indicating its relevance, then the Discounted Cumulative Gain at K is
$DCG@K=\frac{\sum_{k=1}^{K}r_{k}}{\log_{2}(k+1)}$; the Ideal $DCG@K$
($IDCG@K$) is the maximum achievable $DCG@K$ given a set of relevance scores
and finally $NDCG@K=\frac{DCG@K}{IDCG@K}$. The definitions of $AP@K$ and
$NDCG@K$ given above are based on the concept of relevancy of the result,
which is not immediate to state in our scenario. Intuitively, we consider a
retrieved formula $\hat{\varphi}$ semantically similar to a searched formula
$\varphi$ if they behave similarly across an arbitrary number of trajectory,
i.e. if they are (un)satisfied on the same subset of signals. Hence, given a
similarity threshold $\omega\in[0,1]$ and a set of trajectories
$\\{\xi_{j}\\}_{j=1}^{n_{\text{traj}}}\subset\mathcal{T}$, we say that
$\hat{\varphi}$ is a relevant result w.r.t. the query $\varphi$ if
$\frac{|\\{j:s(\varphi,\xi_{j})=s(\hat{\varphi},\xi_{j})\\}|}{n_{\text{traj}}}\geq\omega$.
Moreover, we also consider the kernel similarity of Equation 2 (whose maximum
is $1$) as an indicator of relevancy of the results, which accounts also for
the quantitative semantics of formulae. In our experiments, we used
$n_{\text{traj}}=$10\,000$$ signals sampled from $\mu_{0}$, $\omega=0.9$ and
STL formulae sampled from $\mathcal{F}$ with $N=3$, varying the number
$n_{\text{nodes}}$ of nodes allowed in the syntax tree of the formula. Results
in terms of $AP@K$ and $NDCG@K$ (with $K=5$) are reported in Table 1,
highlighting both the ranking quality of our methodology ($NDCG@K$ always $1$)
and the relevancy of retrieved results (median<EMAIL_ADDRESS>in the worst case),
as well as resilience to the size of the searched formulae (we recall that the
database contains specification with $n_{\text{nodes}}\leq 5$). Moreover, we
also observe that the kernel similarity between the query $\varphi$ and the
top retrieved formula $\hat{\varphi}$ is high (median above $0.84$ in all
tested cases) across all specification dimensions.

###  3.2 Bayesian Optimization in the Semantic Space of Formulae

Algorithm 3 STL Requirement Mining Combining BO and IR $\mathcal{D}_{p}$,
$\mathcal{D}_{n}$, SemanticDB  $\tilde{\varphi}$ k_phis, phis, g_phis
$\leftarrow getInitialCandidates(\mathcal{D}_{p},\mathcal{D}_{n})$
$i\leftarrow 0$ opt$\leftarrow[0,g\\_phis]$ while $i\leq maxiter$ or
$(opt[i]-opt[i-1])<\varepsilon$ do new_phis$\leftarrow$optimizeUCB()
k_phis.append(new_phis)  phis.append(retrieveFromEmbedding(new_phis,
SemanticDB))  g_phis.append(evaluteObjective(new_phis, $\mathcal{D}_{p}$,
$\mathcal{D}_{n}$))  opt.append($\max(g\\_phis)$)  end while
$\tilde{\varphi}\leftarrow$ phis[argmax(g_phis)]

As already mentioned, the goal of the methodology we propose is to learn the
function mapping semantic embeddings of STL formulae (built as for Equation 2)
to values of $G(\varphi)$ of Equation 3, given labeled datasets
$\mathcal{D}_{p}$ and $\mathcal{D}_{n}$ of respectively positive and negative
trajectories. This implies searching in the latent semantic space for the
continuous representation of a formula able to discriminate between signals
belonging to $\mathcal{D}_{p}$ and $\mathcal{D}_{n}$. Being it a non-linear
non-convex optimization problem, we tackle it by means of Bayesian
Optimization (BO), and more specifically of the Gaussian Process Upper-
Confidence Bound (GP-UCB) algorithm [26]. We recall that the typical BO loop
consists in: (i) fit the GP on the available data (pairs
$(k(\varphi),G(\varphi))$ in our case, with $\varphi$ STL formula and
$k(\varphi)$ its STL kernel embedding); (ii) optimize the UCB acquisition
function to find the best candidate (the kernel embedding $k(\hat{\varphi})$
of a formula $\hat{\varphi}$ improving Equation 3); (iii) query the objective
function, i.e. Equation 3, on the new points to add new data pairs to the
training set. In order to obtain an STL formula $\hat{\varphi}$ starting from
its emebedding $k(\hat{\varphi})$ (i.e. between the steps (ii) and (iii) of
the BO loop), we augment the procedure with a retrieval step, querying the
semantic vector database on $k(\hat{\varphi})$. An overview of the overall
methodology is reported in Algorithm 3. In our experiments we use GP with
Matern kernel (starting with an initial batch of only $10$ points) and
optimize the UCB acquisition function using Stochastic Gradient Descent (SGD),
given the high-dimensionality of the latent semantic space (namely $1000$, as
detailed in Section 3.1); this is done in Python exploiting the BoTorch
library [4], which allows to leverage GPU acceleration.

##  4 Experiments

###  4.1 Experimental Setting

We test our proposed algorithm across a number of benchmark datasets found in
the related literature (see Section 5). Unless differently specified, we
perform $5$-fold cross validation and report mean and standard deviation of
the results across the folds. The metrics we use for assessing the quality of
the results are: Misclassification Rate ($MCR$), Precision ($Prec$) and Recall
($Rec$), defined as follows, given positive and negative trajectory datasets
$\mathcal{D}_{p}$ and $\mathcal{D}_{n}$ and the mined STL specification
$\hat{\varphi}$:

$\displaystyle MCR$ $\displaystyle=$
$\displaystyle\frac{|\\{\xi_{j}\in\mathcal{D}_{p}:s(\hat{\varphi},\xi_{j})<0\\}|+|\\{\xi_{j}\in\mathcal{D}_{n}:s(\hat{\varphi},\xi_{j})>0\\}|}{|\mathcal{D}_{p}|+|\mathcal{D}_{n}|}$
$\displaystyle Prec$ $\displaystyle=$
$\displaystyle\frac{|\\{\xi_{j}\in\mathcal{D}_{p}:s(\hat{\varphi},\xi_{j})>0\\}|}{|\\{\xi_{j}\in\mathcal{D}_{p}:s(\hat{\varphi},\xi_{j})>0\\}|+|\\{\xi_{j}\in\mathcal{D}_{n}:s(\hat{\varphi},\xi_{j})>0\\}|}$
$\displaystyle Rec$ $\displaystyle=$
$\displaystyle\frac{|\\{\xi_{j}\in\mathcal{D}_{p}:s(\hat{\varphi},\xi_{j})>0\\}|}{|\\{\xi_{j}\in\mathcal{D}_{p}:s(\hat{\varphi},\xi_{j})>0\\}|+|\\{\xi_{j}\in\mathcal{D}_{p}:s(\hat{\varphi},\xi_{j})<0\\}|}$

In the above Equation, using the jargon of binary classification algorithms,
the true positives are
$TP=\\{\xi_{j}\in\mathcal{D}_{p}:s(\hat{\varphi},\xi_{j})>0\\}$, the true
negatives are $TN=\\{\xi_{j}\in\mathcal{D}_{n}:s(\hat{\varphi},\xi_{j})<0\\}$,
false positives are instead defined as
$FP=\\{\xi_{j}\in\mathcal{D}_{n}:s(\hat{\varphi},\xi_{j})>0\\}$ and finally
false negatives
$FN=\\{\xi_{j}\in\mathcal{D}_{p}:s(\hat{\varphi},\xi_{j})<0\\}$.

The retrieval step of Algorithm 3 requires to normalize trajectories of
$\mathcal{D}_{p}$ and $\mathcal{D}_{n}$ by removing the mean and diving by the
standard deviation of each dimension, in order to align their domain with that
of $\mu_{0}$ (i.e. the measure w.r.t. which semantic embeddings are computed).
At the end of the procedure, they are scaled back to their original range, and
the same inverse normalization is applied to the thresholds of the variables
appearing in the mined specification (since all the formulae in the semantic
vector database have thresholds compatible with $\mu_{0}$ by construction).
Another aspect to take into account when dealing with the retrieval step of
Algorithm 3 is that of time parameters: the sampling rate of $\mu_{0}$ is $1$
per time-step, considering a total of $100$ points, but this might vary across
trajectory distribution. Assuming that the target distribution
$\mathcal{D}_{\text{test}}$ has $n_{\text{test}}$ time points, without loss of
generality we can consider them to be sampled at rate of $1$ per time step, so
that a time interval $[a,b]$ referred to a trajectory sampled from $\mu_{0}$
can be mapped to the corresponding time interval
$[a_{\text{test}},b_{\text{test}}]$ for a trajectory sampled from
$\mathcal{D}_{\text{test}}$ as:
$[a_{\text{test}},b_{\text{test}}]=[\left\lfloor
a\cdot(n_{\text{test}}/100)\right\rfloor,\left\lfloor
a\cdot(n_{\text{test}}/100)\right\rfloor+\left\lceil(b-a)\cdot(n_{\text{test}}/100)\right\rceil]$.
Moreover, before querying the semantic vector database, we restrict the search
space to only those indexes containing formulae having a number of variables
compatible with the dimension of $\mathcal{D}_{\text{test}}$ and, for the sake
of interpretability, we first search on indexes containing specifications with
at most $4$ nodes in the syntactic tree and, if the outcome of the procedure
is not satisfactory, we allow also to search for more complex requirements (as
we will report in the remainder of this Section, this has not been necessary
in any of test cases). All the experiments are implemented in Python,
leveraging GPU acceleration in every sub-procedure thanks to the PyTorch [22]
library, and tested on a AMD EPYC $7542$ machine with $64$ GB of RAM and a
NVidia A$100$ with $80$GB of memory.

###  4.2 Case Studies

(a) Linear System (b) HAR (c) Control of a train (d) Maritime (e) LP$5$

Figure 1: Regular and anomalous trajectories sampled from the testes datasets.

##### Linear System

In [18] the authors generate trajectories using the following dynamical
system: regular trajectories evolve as $\dot{x}=0.03\cdot x+w$ while anomalous
ones as $\dot{x}=-0.03\cdot x+w$, being $w(t)$ a white noise with variance
$0.04$. Following their approach, we generated $100$ positive and $100$
negative examples, sampling $100$ points for each trajectory, a random subset
of the obtained signals is reported in Figure 1(a).

##### Maritime surveillance

The Maritime Surveillance dataset is a common benchmark for STL anomaly
detection, used e.g. in [8, 18, 19]. It consists in $2$-dimensional signals
reporting vessels behaviour, divided in $1000$ regular and $1000$ anomalous
examples, each consisting of $61$ sampling points. Random signals belonging to
the maritime dataset are reported in Figure 1(d).

##### Human Activity Recognition (HAR)

The HAR dataset created in [1] and publicly available at [14] contains
recording of humans performing activities in daily livings. Following the
approach [17], we split the $6$ possible activities in $3$ static postures and
$3$ dynamic activities, and use these as classification labels (for a total of
$36$ and $58$ trajectories belonging to the mentioned classes, respectively).
Signals belonging to this dataset are $1$-dimensional, representing the
standard deviation of the acceleration of the subjects in the $x$-direction at
$50$ sampling points. A subset of positive and negative trajectories belonging
to the HAR dataset is reported in Figure 1(b).

##### Robot Execution Failures in Motion with Part (LP$5$)

The LP$5$ dataset, available at [14], consists of $3$-dimensional signals,
representing force measurements along the $3$ axes, each recorded at $15$ time
steps. As in [17], we consider as positive the $44$ instances showing normal
behaviour and as negative the $26$ exhibiting bottom collisions. Trajectories
belonging to this dataset are shown in Figure 1(e).

##### Cruise Control of Train

Following [18], we test our procedure on the dataset coming from the (noisy)
observations of the velocity of a train, simulated in a cruise control
scenario. Such signals are $1$-dimensional and sampled for $100$ points;
classes contain $100$ trajectories each. A random sample of such trajectories
is shown in Figure 1(c).

Table 2: Best mined formula, mean and std of MCR, Prec and Rec in the test set across the folds for each case study, computational time (in seconds) required to find the reported formula using Algorithm 3. $\hat{\varphi}$ | $MCR$ | $Prec$ | $Rec$ | Time ($s$)
---|---|---|---|---
Linear | $F_{[70,\infty]}(x_{0}\leq 1.16)$ | $0.00\pm 0.00$ | $1.00\pm 0.00$ | $1.00\pm 1.00$ | $7.32\pm 1.45$
Maritime | $(x_{1}\geq 23.19)U(x_{0}\leq 32.56)$ | $0.07\pm 0.01$ | $0.97\pm 0.01$ | $1.00\pm 0.00$ | $12.4\pm 2.63$
HAR | $G(x_{0}\geq-0.75)$ | $0.00\pm 0.00$ | $1.00\pm 0.00$ | $1.00\pm 0.00$ | $8.32\pm 0.25$
LP$5$ | $G_{[4,\infty]}(x_{2}\geq 0.084)$ | $0.01\pm 0.01$ | $0.98\pm 0.01$ | $1.00\pm 0.00$ | $11.6\pm 1.44$
Train | $G_{[0,36]}(x_{0}\leq 37)$ | $0.03\pm 0.01$ | $0.93\pm 0.02$ | $1.00\pm 0.00$ | $8.21\pm 1.23$

####  4.2.1 Discussion

Figure 2: Signals sampled from $\mu_{0}$.

Results of the retrieval-augmented STL requirement mining procedure are shown
in Table 2, where we report the best formula we infer (denoted as
$\hat{\varphi}$), along with performance metrics and computational time for
each case study. In this context, with best formula we denote the mined
formula having the lowest $MCR$ in the test-set, using the lower number of
nodes as tiebreaker. Notably, we can observe that all mined requirements are
interpretable, having a syntax tree of at most $3$ nodes, in line with more
recent related works (namely [17, 18, 19]), and act as almost-perfect
discriminators among regular and anomalous trajectories (see $MCR$ in Table
2). Moreover, we can observe that the mined specifications never produce false
positives (reporting a recall of $1$ in every test case). It is worth noting
that even if the LP$5$ benchmark is $3$-dimensional, the mining procedure
outputs a formula in which only a variable appears (in line with [17]). We
find it interesting to remark the difference between trajectories sampled from
the measure $\mu_{0}$ used to compute STL embeddings (a random sample of
signals belonging to this distribution is shown in Figure 2) and those of the
test benchmarks (reported in Figure 1), to underline the applicability of our
methodology across various trajectories distribution, without the need of
recomputing the semantic vector database on each data distribution. Finally,
another relevant difference between our approach and related works is that we
fully exploit GPU acceleration in every step of the mining procedure,
resulting in lower computational times. Indeed as reported in Table 2, our
method is able to learn a classifier in less than $15s$ for all the tested
cases, while [18] reports a computational time of $39.05s,23.17s,32.32s$ for
the Linear, Maritime and Train datasets, respectively; [17] register $22.07s$
and $22.83s$ for mining a discriminating formula in the HAR and LP$5$
banchmarks; finally [19] takes $73s$ to classify signals belonging to the
Maritime dataset.

##  5 Related Work

##### Mining STL specifications from data

has seen a surge of interest in the last few years, as reported in [7]. At a
high level, following [7], contributions in this field can be characterized
by: template-based vs template-free procedures (in the first, the mining
consists in finding the optimal parameters for a given syntactic skeleton,
while in the second both the structure and the parameters of the specification
have to be inferred), supervised vs unsupervised approaches (following the
availability of labeled trajectory datasets) and online vs offline methodology
(the first concerns updating the inferred requirements every time new data are
available, the latter considers instead a static set of data). In this
scenario our algorithm can be defined as a template-free, supervised, offline
approach (besides, it can be cast as an unsupervised learning problem by just
changing the objective function of Equation 3, keeping intact the structure of
Algorithm 3). Most of the related works with these characteristics decompose
the task in two steps, i.e. as bi-level optimization problem: they first learn
the structure of the specification from data, and then instantiate a concrete
formula using parameter inference methods. For the first step, a popular
extension of STL, called Parametric Signal Temporal Logic (PSTL) has been
devised [2], in which parameters replace both threshold constants in numerical
predicates and time bounds in temporal operators. In [27] STL specifications
are mined using a grid-based discretization of both time and value domain,
used to guide a clusterization of the input signals, for assigning to each
cluster a descriptive STL formula and obtain the final requirement as
disjunction of cluster formulae. In [8] a map between a fragment of PSTL and
binary decision trees is established, so that the requirement mining procedure
is cast into a classification problem tackled by decision tree, and then
translation of the tree to a formula is done to infer the final specification.
In [18] STL requirements are mined using systematic enumeration of PSTL
templates (in increasing order w.r.t. the syntactic size) and accuracy-driven
instantiation of each template (towards optimizing logical classification),
guaranteeing that the inferred formula is the smallest one reaching a given
accuracy threshold. Another line of works [3, 19, 25] use instead Genetic
Algorithms (GA) to learn the PSTL structure of the specification, and find
optimal parameters by Gaussian Process Upper Confidence Bound (GP-UCB);
notably [19] do not put any restrictions on the fragment of PSTL used. Our
work differentiate from the cited papers in that it learns simultaneously both
the structure and the parameters of the STL requirement, by performing
optimization in a continuous space representing the semantics of formulae.

##### Data Mining for Time Series Data

Typically, anomaly detection and classification in time series data is done by
either autoregressive statistical models, distance and density-based
clustering methods, or prediction and reconstruction-based deep neural
networks [21, 11]. Such methods often require to individuate a latent
representation of trajectories or to hand-craft features for the analysis,
leading to a black-box solution to the problem, which do not offer any insight
on the system at hand. Notably, there are works focusing on building time
series classifiers which are interpretable, e.g. [10] that uses Symbolic
Fourier Approximation to infer a symbolic representation of the input before
learning a linear classifier, and [20] that uses saliency maps to characterize
classes. Using STL as formal language for the learned class discriminator
allows to describe emergent properties of a system in a rich yet interpretable
way, leading to knowledge discovery about the process under analysis [7].

##### Information Retrieval and Data Mining,

as well as their interaction, witnessed a growth of interest recently [15]. In
particular, the adoption of the so-called dense vector databases (i.e.
specialized storage systems designed for efficient management of dense vectors
and supporting advanced similarity search) and retrieval methods has been
proven to be effective in injecting factual information into Machine Learning
(ML) techniques [28]. Using information retrieval (IR) systems to support ML
models access abstractly-represented knowledge stores helps them in generalize
beyond training data, without the need of increasing their size, and notably
helps for the interpretability and explainability of the outcomes, being
inference grounded on retrieved information, often stored in a human-readable
format [28].

##  6 Conclusion

In this work we propose a methodology for mining STL requirements from time-
series data, in a purely data-driven setting. In particular, we focus on
binary classification of trajectories, such as those coming from observations
of Cyber-Physical Systems. The core idea of our algorithm lies in the
interaction between Bayesian Optimization and Information Retrieval
techniques, by searching from the discriminating STL property directly in a
continuous space representing its robust semantics. Doing so, we learn
simultaneously the structure and the parameters of the specification, in such
a way that it is possible to control its maximum allowed size, enhancing the
interpretability of the inferred formula, to foster knowledge discovery about
the system. Results of experiments on several benchmarks show that we achieve
state-of-art accuracy, drastically reducing computational time w.r.t. related
works, thanks to the exploitation of GPU acceleration in all steps of our
methodology. It is also worth mentioning that by changing the optimization
objective, our framework can be deployed in an unsupervised classification
scenario, as well as in a requirement mining context in which only positive
examples are available, extensions that we plan to implement as future work.

## References

  * [1] Anguita, D., Ghio, A., Oneto, L., Parra, X., Reyes-Ortiz, J.L.: In: 21st European Symposium on Artificial Neural Networks, ESANN 2013 (2013)
  * [2] Asarin, E., Donzé, A., Maler, O., Nickovic, D.: Parametric identification of temporal properties. In: Runtime Verification - Second International Conference, RV 2011, Revised Selected Papers. Lecture Notes in Computer Science, vol. 7186, pp. 147–160. Springer (2011)
  * [3] Aydin, S.K., Gol, E.A.: Synthesis of monitoring rules with STL. J. Circuits Syst. Comput. 29(11), 2050177:1–2050177:26 (2020)
  * [4] Balandat, M., Karrer, B., Jiang, D.R., Daulton, S., Letham, B., Wilson, A.G., Bakshy, E.: BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization. In: Advances in Neural Information Processing Systems 33 (2020)
  * [5] Bartocci, E., Bortolussi, L., Nenzi, L., Sanguinetti, G.: System design of stochastic models using robustness of temporal properties. Theor. Comput. Sci. 587, 3–25 (2015)
  * [6] Bartocci, E., Deshmukh, J.V., Donzé, A., Fainekos, G., Maler, O., Nickovic, D., Sankaranarayanan, S.: Specification-based monitoring of cyber-physical systems: A survey on theory, tools and applications. In: Lectures on Runtime Verification - Introductory and Advanced Topics, Lecture Notes in Computer Science, vol. 10457, pp. 135–175. Springer (2018)
  * [7] Bartocci, E., Mateis, C., Nesterini, E., Nickovic, D.: Survey on mining signal temporal logic specifications. Inf. Comput. 289(Part), 104957 (2022)
  * [8] Bombara, G., Vasile, C.I., Penedo, F., Yasuoka, H., Belta, C.: A decision tree approach to data classification using signal temporal logic. In: Proceedings of the 19th International Conference on Hybrid Systems: Computation and Control, HSCC 2016. pp. 1–10. ACM (2016)
  * [9] Bortolussi, L., Gallo, G.M., Kretínský, J., Nenzi, L.: Learning model checking and the kernel trick for signal temporal logic on stochastic processes. In: Tools and Algorithms for the Construction and Analysis of Systems - 28th International Conference, TACAS 2022,Proceedings, Part I (2022)
  * [10] Early, J., Cheung, G.K.C., Cutajar, K., Xie, H., Kandola, J., Twomey, N.: Inherently interpretable time series classification via multiple instance learning. CoRR abs/2311.10049 (2023)
  * [11] Faouzi, J.: Time Series Classification: A review of Algorithms and Implementations. In: Machine Learning (Emerging Trends and Applications). Proud Pen (2022)
  * [12] Han, Y., Liu, C., Wang, P.: A comprehensive survey on vector database: Storage and retrieval technique, challenge (2023)
  * [13] Johnson, J., Douze, M., Jégou, H.: Billion-scale similarity search with GPUs. IEEE Transactions on Big Data 7(3), 535–547 (2019)
  * [14] Kelly, M., Longjohn, R., Nottingham, K.: The uci machine learning repository, https://archive.ics.uci.edu
  * [15] Liu, J., Kong, X., Zhou, X., Wang, L., Zhang, D., Lee, I., Xu, B., Xia, F.: Data mining and information retrieval in the 21st century: A bibliographic review. Comput. Sci. Rev. 34 (2019)
  * [16] Maler, O., Nickovic, D.: Monitoring temporal properties of continuous signals. In: Proc. of FORMATS and FTRTFT, LNCS, vol. 3253, Springer (2004), pp. 152-166. Lecture Notes in Computer Science, vol. 3253, pp. 152–166. Springer (2004)
  * [17] Mohammadinejad, S., Deshmukh, J.V., Puranic, A.G., Vazquez-Chanlatte, M., Donzé, A.: Interpretable classification of time-series data using efficient enumerative techniques (2019)
  * [18] Mohammadinejad, S., Deshmukh, J.V., Puranic, A.G., Vazquez-Chanlatte, M., Donzé, A.: Interpretable classification of time-series data using efficient enumerative techniques. In: HSCC ’20: 23rd ACM International Conference on Hybrid Systems: Computation and Control. pp. 9:1–9:10. ACM (2020)
  * [19] Nenzi, L., Silvetti, S., Bartocci, E., Bortolussi, L.: A robust genetic algorithm for learning temporal specifications from data. In: Quantitative Evaluation of Systems - 15th International Conference, Beijing, China, Proceedings. Lecture Notes in Computer Science, vol. 11024, pp. 323–338. Springer (2018)
  * [20] Nguyen, T.L., Gsponer, S., Ilie, I., O’Reilly, M., Ifrim, G.: Interpretable time series classification using linear models and multi-resolution multi-domain symbolic representations. CoRR abs/2006.01667 (2020)
  * [21] Oswal, S., Shinde, S., Vijayalakshmi, M.: A survey of statistical, machine learning, and deep learning-based anomaly detection techniques for time series. In: Advanced Computing. pp. 221–234. Springer Nature Switzerland, Cham (2023)
  * [22] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A.: Automatic differentiation in pytorch. In: NIPS-W (2017)
  * [23] Pnueli, A.: The temporal logic of programs. In: 18th Annual Symposium on Foundations of Computer Science (sfcs 1977). pp. 46–57 (1977)
  * [24] Rasmussen, C.E., Williams, C.K.I.: Gaussian processes for machine learning. Adaptive computation and machine learning, MIT Press (2006), https://www.worldcat.org/oclc/61285753
  * [25] Saglam, I., Gol, E.A.: Cause mining and controller synthesis with STL. In: 58th IEEE Conference on Decision and Control, CDC 2019. pp. 4589–4594. IEEE (2019)
  * [26] Srinivas, N., Krause, A., Kakade, S.M., Seeger, M.W.: Information-theoretic regret bounds for gaussian process optimization in the bandit setting. IEEE Trans. Inf. Theory 58(5), 3250–3265 (2012)
  * [27] Vaidyanathan, P., Ivison, R., Bombara, G., DeLateur, N.A., Weiss, R., Densmore, D., Belta, C.: Grid-based temporal logic inference. In: 56th IEEE Annual Conference on Decision and Control, CDC 2017. pp. 5354–5359. IEEE (2017)
  * [28] Zamani, H., Diaz, F., Dehghani, M., Metzler, D., Bendersky, M.: Retrieval-enhanced machine learning. In: SIGIR ’22: The 45th International ACM SIGIR Conference on Research and Development in Information Retrieval. pp. 2875–2886. ACM (2022)