# Reconstructed Student-Teacher and Discriminative Networks
for Anomaly Detection

Shinji Yamada1, Satoshi Kamiya2 and Kazuhiro Hotta3 1Meijo University, 1-501
Shiogamaguchi, Tempaku-ku, Nagoya 468-8502, Japan,
<EMAIL_ADDRESS>University, 1-501 Shiogamaguchi,
Tempaku-ku, Nagoya 468-8502, Japan<EMAIL_ADDRESS>University, 1-501 Shiogamaguchi, Tempaku-ku, Nagoya 468-8502, Japan,
<EMAIL_ADDRESS>
###### Abstract

Anomaly detection is an important problem in computer vision; however, the
scarcity of anomalous samples makes this task difficult. Thus, recent anomaly
detection methods have used only ”normal images” with no abnormal areas for
training. In this work, a powerful anomaly detection method is proposed based
on student-teacher feature pyramid matching (STPM), which consists of a
student and teacher network. Generative models are another approach to anomaly
detection. They reconstruct normal images from an input and compute the
difference between the predicted normal and the input. Unfortunately, STPM
does not have the ability to generate normal images. To improve the accuracy
of STPM, this work uses a student network, as in generative models, to
reconstruct normal features. This improves the accuracy; however, the anomaly
maps for normal images are not clean because STPM does not use anomaly images
for training, which decreases the accuracy of the image-level anomaly
detection. To further improve accuracy, a discriminative network trained with
pseudo-anomalies from anomaly maps is used in our method, which consists of
two pairs of student-teacher networks and a discriminative network. The method
displayed high accuracy on the MVTec anomaly detection dataset.

##  I INTRODUCTION

Anomaly detection is the identification of samples that deviate from normal
patterns in the data. In recent years, convolutional neural networks (CNNs)
have been applied to anomaly detection for the visual inspection of industrial
products[2, 3, 4], surveillance of locate intruders[5, 6, 7], and pathological
diagnosis of medical images[8, 9, 10, 11, 12]. In the past, image-level
anomaly detection methods that classify each image as an anomaly have been
studied, but image-level anomaly detection methods are not suitable for the
inspection of products that require anomaly detection at the pixel level.
Famous image-level anomaly detection methods rely on generative models to
accurately reconstruct normal images. The generative model, based on the
generative adversarial network (GAN)[11] and autoencoders [13], can be trained
using only normal images which allow the reconstruction of normal regions
although anomalous regions cannot be reconstructed because the model is
trained using only normal images. Anomaly detection using generative models
identifies anomalous regions by visualizing regions with poor reconstruction
capabilities. However, it is difficult for these methods to identify anomalous
regions unless the generated normal regions can be reconstructed with high
accuracy[15, 13].

Figure 1: The problem of STPM-based methods. The anomaly map has large values
for both abnormal and normal regions.

The use of pretrained models has been proposed as a new approach for anomaly
detection, such as STPM[14] which uses a student-teacher network. In STPM, the
teacher network is the pretrained ResNet18, and the student network is the
untrained ResNet18. Only normal images are used for STPM training. The student
network learns, thereby making feature maps similar to those in the teacher
network. Because training is performed only on normal images, the student
network can only output the features of normal regions. In contrast, the
teacher network is a model pre-trained on ImageNet, so it can represent
features of abnormal regions well. The difference in the feature
representation between the student and teacher networks is the anomaly region.

In student ResNet18, STPM detects anomalous regions using anomaly maps at
three different resolutions. However, STPMs have several limitations. STPM
uses the difference between feature maps at three different resolutions in
student and teacher networks. In the validation phase, the final anomaly map
is calculated by multiplying the three maps. When all three anomaly maps
detect the same abnormal area, highly accurate abnormality detection is
possible. However, if one of the three maps fails to detect an abnormal
region, the anomaly region cannot be detected because the three maps are
multiplied.

Therefore, each of the three anomaly maps must be improved. Observing that
STPM does not have the ability to reconstruct normal images similar to
generative models which are effective in anomaly detection, we introduced a
new student network that has the ability to reconstruct normal features.
Conventional STPM uses ResNet18 as the teacher network. If the same teacher
network is used for a new student, it is possible that similar anomaly maps
can be obtained. To obtain anomaly maps from different viewpoints, the pre-
trained ResNet50 was used as the teacher network for the student network for
reconstruction. However, because the architecture of the student network is
ResNet18, knowledge distillation for student-teacher pairs with different
structures is more challenging than learning student-teacher pairs as in STPM
with the same structure. Therefore, in the proposed method, to proceed with
the learning, an attention mechanism is used with the structure to propagate
some features of the teacher network to the student network. Because only
normal images are used for training in the proposed method, the attention
mechanism reconstructs normal regions with higher accuracy. Knowledge
distillation is ensured by passing hints through the attention mechanism.

Figure 2: STPM and its problems. STPM detects anomalies by multiplying the
anomaly maps at three different resolutions. However, if one of the three
anomaly maps fails to detect an anomaly location, the anomaly detection fails.
The bottom images show examples of failures.

Although the use of a student-teacher pair for reconstruction improves the
accuracy of anomaly detection, a discriminative network is added to further
improve the accuracy of anomaly detection. Because anomaly detection with
knowledge distillation uses only normal images for training, it cannot
discriminate whether a product is an anomaly. Thus, there are cases in which
anomaly maps have large values for the normal regions, as shown in Figure 1.
Anomaly maps can detect abnormal areas, but they also detect normal areas as
anomalies. This decreases the accuracy of anomaly detection.

Therefore, the anomaly maps obtained by our STPM, with two pairs of student-
teacher networks, are input into a discriminative network, to reconsider the
anomaly map when pseudo-abnormal images are also fed into our STPM. The
discriminative network thus learns to produce a more accurate anomaly map. The
accuracy of anomaly detection is improved by multiplying the anomaly map
obtained from the discriminant network by that obtained from the STPM with two
pairs of networks.

The proposed method was evaluated on an MVTec anomaly detection dataset[2].
The pixel-level AUC of the proposed method outperformed that of the original
STPM in many categories. The image-level AUC of the proposed method was also
better than that of the original STPM. In addition, the effectiveness of
different teacher networks and attention mechanisms was demonstrated in
ablation studies.

The rest of this paper is organized as follows. In Section II, related works
are described. In Section III, the proposed method is explained in detail. In
Section IV, the experimental results are presented. Finally, the study is
concluded in Section V.

Figure 3: Overview of STPM with two pairs of student-teacher networks.

##  II Related Works

###  II-A Anomaly Detection

There are two types of anomaly detection: image-level and pixel-level. The
goal of image-level anomaly detection is to classify anomalous samples
correctly. There are three approaches for image-level anomaly detection:
generative models, the distribution of data such as feature spaces, and
classification. The generative model detects anomalies based on the degree of
reconstruction loss[16, 17, 18]. The distribution-based method considers
samples that deviate from the normal data distribution as abnormal. When
generating a probability distribution of only normal products, the probability
density of abnormal products is low[19, 20], allowing the classification of
abnormal data. A method based on classification is a method for anomaly
detection that combines geometric transformation and classification[21]. It
detects anomalous samples based on the idea that classification accuracy for
unknown anomalous data is poor. However, these methods are unsuitable for
anomaly detection at the pixel level.

Pixel-level anomaly detection is a method used for the detection of anomalies
at each pixel. The difficulty of detecting anomalies at the pixel level is
higher than that of detecting anomalies at the image level because the number
of inspection targets increases. The mainstream method for detecting anomalies
at the pixel level is based on generative models[13, 10, 11, 22], such as GAN
and autoencoders. In recent years, methods that combine GAN and autoencoders
have been studied[23, 24, 25]. The generative model must reconstruct normal
images with high accuracy; otherwise, the accuracy of anomaly detection
decreases.

SPADE [26] was proposed as a new approach for pixel-level anomaly detection.
SPADE compares the features of normal and abnormal images based on a pre-
trained model, and k-means clustering was used to detect abnormalities. In
addition to the pre-trained model, uninformed students [1] using knowledge
distillation with a student-teacher structure have also been proposed.
Uninformed students train the student network using only normal data. Anomaly
maps are then calculated from two perspectives: (i) the difference between the
outputs of student and teacher networks, (ii) the uncertainty of multiple
student networks. STPM[14] was proposed as a method for developing both SPADE
and uninformed students. In STPM, because the training is performed only on
normal data, the student network represents only the features of the normal
region. In contrast, the teacher network is pre-trained on ImageNet, so it can
represent the features of anomalous regions. The difference in the feature
representation between the student and teacher networks is the anomaly region.
STPM uses knowledge distillation for feature maps at three different
resolutions and outputs three anomaly maps. The final anomaly map is computed
by multiplying the three anomaly maps. However, if one of the three anomaly
maps is unable to detect an anomalous region, the final anomaly map cannot
detect the anomalous region, as shown in Figure 2. To address this problem,
each anomaly map is improved by using a new student network for
reconstruction.

###  II-B Attention Mechanism

Various attention mechanisms have been proposed for use in image recognition.
The residual attention network [27] solved the vanishing gradient problem
using a structure similar to that of a residual block. Squeeze-and-excitation
network (SENet)[28] introduced an attention mechanism that emphasizes
important channels in feature maps. A transformer [29] was proposed for
language translation using only the attention mechanism. Several image
recognition methods using self-attention, a type of transformer, have also
been proposed[30, 31, 32]. An attention branch network[33] proposes an
attention map for classification by aggregating multiple feature maps. An
attention map can be used to visualize the basis of decisions.

In our method, the purpose of the attention mechanism is to leak the features
of the teacher network to the student network to effectively reconstruct the
features of normal data. Because anomalies must be detected at the pixel
level, an attention mechanism that can emphasize and suppress pixels is more
suitable than an attention mechanism that emphasizes channels such as SENet.
If almost all features in the teacher network are leaked into the student
network, there will be no difference between the student and teacher networks.
Therefore, by aggregating the features in the teacher network into one
channel, pixels can be emphasized and suppressed without providing the entire
information. Therefore, we used an attention map generated from the teacher
network. Because only normal data are used for training the student network,
the attention mechanism is used to reconstruct the normal regions.

##  III Proposed Method

The details of the proposed method are described in this section. The proposed
method makes four contributions: the student network for reconstruction,
attention from the teacher network to the student network, different teacher
networks from the original STPM, and the discriminative network.

In Section III.A, the original STPM and its problems, in Section III.B, the
student network for reconstruction and its teacher network, in Section III.C,
the attention mechanism from the teacher to the student network, and in
Section III.D, a discriminative network that revisits anomaly maps obtained by
STPM with two pairs of student-teacher networks are described.

###  III-A Problem of STPM

STPM is an effective anomaly detection method. STPM multiplies the three
anomaly maps at different resolutions. The accuracy of the anomaly detection
was improved by multiplying the three anomaly maps. However, if one of the
three anomaly maps is unable to detect the anomalous region, the anomalous
region cannot be detected after multiplying the three maps. Figure 2 shows an
example in which STPM failed to detect anomalous regions because of such a
problem. Multiplying the three anomaly maps improved the accuracy of anomaly
detection with respect to the detection accuracy at each resolution, but there
were cases of failure after multiplication. To address this problem, it was
necessary to improve each anomaly map at three different resolutions. As shown
in the network structure of STPM in Figure 2, STPM does not use feature
reconstruction such as generative models. Adding reconstruction abilities,
such as generative models, to new student networks was considered, to improve
the accuracy of anomaly detection. The proposed network detects anomalies
using two pairs of student-teacher networks.

Figure 4: Attention mechanism in the proposed method. The student network can
emphasize the important pixels for reconstruction. Figure 5: Visualization of
attention map.  Figure 6: Overview of the proposed network. The discriminant
network uses pseudo-anomalies to reconsider the anomaly map obtained from STPM
with two pairs of student-teacher networks. In the test phase, the anomaly map
from the STPM is multiplied by that from the discriminant network.

###  III-B Student Network for Reconstruction

In this work, a new student network is introduced to reconstruct normal
features. There are three main contributions in the new student network. The
first is the reconstruction of normal features. Second, a teacher network
different from that of the original STPM is used for the new student network.
The third is the attention mechanism from the teacher to student network, as
explained in the next section.

Figure 3 presents an overview of the proposed method. The input for the new
student network is the feature map at the lowest resolution in the pre-trained
ResNet18, whereby the student network reconstructs the features of normal
data. The structure of the new student network is similar to that of ResNet18.
Therefore, if the same teacher network is used for student2 as in the original
STPM, the same mistakes as in the original STPM will be made. Thus, ResNet50
pre-trained on ImageNet was used as the teacher network for student2 because
ResNet50 has features and structures that are different from ResNet18, which
is teacher1. With a perspective different from student1, such features will
improve anomaly detection accuracy.

The equation for knowledge distillation is the same as that for conventional
STPM. A normal image is defined as $I_{k}\in\mathbb{R}^{w\times h\times c}$,
where $h$ is the height, $w$ is the width, and $c$ is the number of channels.
The features in the teacher and student networks are normalized along the
channel dimensions as,

$\displaystyle\hat{F_{t}^{l}}{(I_{k})}_{ij}=\frac{F^{l}_{t}(I_{k})_{ij}}{{||F^{l}_{t}(I_{k})_{ij}||}_{\ell_{2}}^{2}},\quad\hat{F_{s}^{l}}{(I_{k})}_{ij}=\frac{F^{l}_{s}(I_{k})_{ij}}{{||F^{l}_{s}(I_{k})_{ij}||}_{\ell_{2}}^{2}}$
(1)

where $F_{t}^{l}$ and $F_{s}^{l}$ represent the features of the teacher and
student networks, respectively; $l$ represents the resolution of the feature
map; and $(i,j)$ is the position. The student network learns to obtain the
same normalized output. The loss function is defined as follows:

$\displaystyle{\ell}^{l}(I_{k})_{ij}=\frac{1}{2}{||\hat{F_{t}^{l}}{(I_{k})}_{ij}-\hat{F_{s}^{l}}{(I_{k})}_{ij}||}_{\ell_{2}}^{2}$
(2)

Equation 2 represents the pixel loss. For the entire image, the loss function
is defined as,

$\displaystyle{\ell}^{l}(I_{k})=\frac{1}{{w_{l}}{h_{l}}}\sum_{i=1}^{w_{l}}\sum_{j=1}^{h_{l}}{\ell^{l}}(I_{k})_{ij}$
(3)

The same loss function is used for feature maps at three resolutions in the
two-student networks. The numbers in Fig. 3 indicate the locations where the
loss function was applied. The final loss function is the sum of the losses,
defined as,

$\displaystyle{\ell}(I_{k})=\sum_{l=1}^{L}{{\ell}^{l}(I_{k})}$ (4)

where L is the number of knowledge distillations. In this study, L = 6 because
the loss was calculated at six locations, as shown in Figure 3. In the test
phase, six anomaly maps $\\{\Omega^{1}(J),...,\Omega^{6}(J)\\}$ were obtained
using student1 and student2. As shown in Figure 3, the anomaly maps at the
same resolution are summed, while the anomaly maps at different resolutions
are multiplied to generate the final anomaly map. When multiplying the three
anomaly maps, the size of each anomaly map is made the same using bilinear
interpolation. Anomaly detection is performed using the final anomaly map.

$\begin{split}\Omega(J)={\\{\Omega^{1}(J)+\Omega^{6}(J)\\}}{\odot}{\\{\Omega^{2}(J)+\Omega^{5}(J)\\}}\\\
{\odot}{\\{\Omega^{3}(J)+\Omega^{4}(J)\\}}\end{split}$ (5)

where $J\in\mathbb{R}^{w\times h\times c}$ denotes the test image. ${\odot}$
is the element-wise multiplication operation. The input for student2 is the
feature maps with a size of 1/32 of the input image, and student2 reconstructs
the feature maps up to 1/4 of the size of the input image. Anomaly maps from
$\Omega^{1}$ to $\Omega^{3}$ in Figure 3 were obtained from the original STPM.
Anomaly maps from $\Omega^{4}$ to $\Omega^{6}$ obtained from the
reconstruction network are also shown. Anomaly map $\Omega^{4}$ has the lowest
resolution. When the anomaly maps of the original STPM are compared with those
of the reconstructed network, the accuracy of the reconstructed network is
better than that of the original STPM. In particular, when the anomaly maps
$\Omega^{1}$ are compared with $\Omega^{6}$ at the highest resolution, it is
observed that the anomaly map $\Omega^{6}$ of the reconstructed network was
able to detect anomalies that were not detected by the original STPM.

###  III-C Attention Mechanism for feature propagation

The third contribution of the proposed method is the attention mechanism from
the teacher network to a new student network for reconstruction. The teacher
network consisted of 50 layers, whereas the student network consisted of 18
layers. Knowledge distillation between networks with different structures is
more difficult to learn than those between networks with the same
structure[34]. It is likely that learning will not proceed as intended. Thus,
an attention mechanism was used to ensure that learning proceeded well. Figure
4 shows the structure of the attention mechanism. The attention mechanism
receives the feature maps in the teacher network and aggregates them into one
channel. The attention mechanism was also used to teach hints to the student
network to facilitate learning. Because only normal data are used for
training, with the attention mechanism, the features in the student network
for the normal region become more similar to the features in the teacher
network. If most of the features in the teacher network are leaked into the
student network, there will be no difference between the student and teacher
networks. Thus, an attention map with only one channel was used to limit the
amount of information provided. Considering an attention mechanism that
emphasizes the channel would not be effective in detecting anomalies because
it would not be able to learn important locations for pixel-level anomaly
detection. Thus, considering a single attention map would be better because it
would emphasize important pixels. The computational graph between the
attention mechanism and the teacher network is not connected; therefore, the
teacher network is not updated. Because only normal products are used as
training data, the attention mechanism learns to successfully reconstruct the
normal regions.

Figure 5 shows examples of attention maps. The attention maps at the highest
resolution, which are feature map 6, are displayed in Figure 3. Figure 5 shows
the attention maps of the normal and abnormal data with three categories along
with the ground truth labels corresponding to the abnormal regions. Attention
maps differ between the normal and abnormal areas. The network was trained
using only the normal data. Therefore, it behaves differently in normal
regions when the input is an abnormal region that has never been seen during
training. The attention mechanism was also introduced into the reconstruction
network in which knowledge distillation is performed. The usage of the
attention mechanism encourages the student network to represent the features
in normal regions.

###  III-D Discriminative Network

In this study, the anomaly detection accuracy of the conventional STPM was
improved by adding a student network for reconstruction. The STPM with two
pairs of student-teacher networks, uses only normal images for training. An
anomaly map is calculated only by taking the difference in features between
the student-teacher networks, and there are cases where the anomaly map has
large values for normal regions. To address this problem, a discriminative
network with U-Net was used. In an ideal anomaly map, the values for the
normal regions are zero and only the abnormal regions have large values. Thus,
a discriminative network was trained to learn pseudo-anomalies to make the
anomaly map closer to the ideal anomaly map.

Figure 6 shows an overview of the proposed method with a discriminative
network. Pseudo-anomalies are input into the trained STPM with two pairs of
student-teacher networks and three anomaly maps are obtained. The anomaly maps
at three resolutions, $\\{\Omega^{1}(I_{a})+\Omega^{6}(I_{a})\\}$,
$\\{\Omega^{2}(I_{a})+\Omega^{5}(I_{a})\\}$, and
$\\{\Omega^{3}(I_{a})+\Omega^{4}(I_{a})\\}$, are fed into the discriminant
network. Note that $I_{a}\in\mathbb{R}^{w\times h\times c}$ is a pseudo-
anomalous image. The discriminative network identifies pseudo-anomalies from
the input anomaly maps. The focal loss used in training the difficult samples
is defined as,

$\displaystyle
L_{seg}(P,T)=\frac{1}{{h}{w}}\sum_{i=1}^{h}\sum_{j=1}^{w}{-T_{ij}(1-{P}_{ij})^{\gamma}\log{{P}_{ij}}}$
(6)

where $P$ is the anomaly map for the pseudo-anomaly predicted by the
discriminative network; $T$ is the ground truth; and $h$ and $w$ denote the
height and width of the input image, respectively. $\gamma$ was set to two.
Because the discriminative network learns to discriminate whether the regions
are anomalous, anomaly values for normal regions are closer to zero, and
anomaly regions have values closer to one. The anomaly map from our STPM was
combined with that of the discriminative network to achieve more accurate
anomaly detection. The final anomaly map was calculated as follows:

$\begin{split}\Omega(J)={\\{\Omega^{1}(J)+\Omega^{6}(J)\\}}{\odot}{\\{\Omega^{2}(J)+\Omega^{5}(J)\\}}\\\
{\odot}{\\{\Omega^{3}(J)+\Omega^{4}(J)\\}}{\odot}{\\{\Omega^{7}(J)\\}}\end{split}$
(7)

where $\Omega^{7}$ is the anomaly map of the discriminant network and $J$ is a
test image. ${\odot}$ is the element-wise multiplication operation.

The pseudo-anomalies proposed in DRAEM[35] were used to train the
discriminative network. DRAEM can generate various types of pseudo-anomalies
using Perlin noise[36] as well as parameters that determine the concentration
of anomalous regions. In this study, DRAEM pseudo-anomalies were used to
prevent overfitting of pseudo-anomalies.

##  IV Experiments

In this section, the experimental results are presented. In Section IV.A, the
datasets and experimental conditions are described. In Section IV.B, the
experimental results for the MVTec dataset are presented. In Section IV.C, the
results of the ablation study are presented.

###  IV-A Datasets

The MVTec anomaly datasets[2] was used for the evaluation. The MVTecAD
includes 15 categories of industrial products. The training data consisted of
only the normal products. The image size was different for each category, and
all images were resized to 256 × 256 pixels. The STPM was trained with two
pairs of student-teacher networks and separately with the discriminative
network. The numbers of epochs for training the STPM and discriminative
network were 100 and 300, respectively. The SGD optimizer was used with a
momentum of 0.9. The learning rate was set to 0.4, batch size to 32, and
weight decay to $1\times 10^{-4}$. The optimization method was Adam in the
discriminative network, and the learning rate was set to 0.0001.

The pixel-level and image-level areas under the ROC curve (AUC) were used as
the evaluation measures. Per-region-overlap (PRO)[1] was also used to evaluate
the anomaly maps.

Figure 7: Examples of anomaly map.

###  IV-B Results on Anomaly detection

First, pixel-level AUC and PRO metrics are described. The pixel-level AUC and
PRO metrics indicate whether the anomaly map correctly detects the anomalous
area. Table I displays the accuracy. Our (Student-Tescher) method in Table I
shows the results using only the network in Figure 3. Our (Discriminative)
method shows the accuracy when only the output $\Omega^{7}$ of the
discriminative network is used. It is observed from the table that the
reconstruction network when combined with the conventional STPM improves the
accuracy, and the discriminative network alone outperforms the conventional
STPM. The anomaly detection accuracy of the proposed method using both the
STPM with two pairs of student-teacher networks and the discriminative
network, as shown in Figure 6, provided further improvement. The pixel-level
AUC and PRO scores of the proposed method outperformed the accuracies of the
DRAEM using pseudo-anomalies and Padim using a pre-trained model. The proposed
method can accurately detect anomalous regions because the reconstruction
network improves the conventional STPM and the discriminative network assists
the STPM using the reconstruction network.

Next, the results of image-level AUC scores are presented. The image-level AUC
is a score that indicates whether the input image can be correctly classified
as abnormal. Table II lists the image-level AUCs. The accuracy of anomaly
detection using the proposed methods outperformed conventional STPM. In the
Table, the anomaly classification accuracies of Padim, STPM, and our (student-
teacher) method which learn only normal products are lower than those of DRAEM
and our (discriminative) method trained using pseudo-anomalies. The accuracy
of anomaly classification can be improved using pseudo-anomalies. The proposed
method achieved the best accuracy by combining the reconstruction and
discriminative networks.

The qualitative results are shown in Figure 7. The anomaly maps of the
proposed method displayed in Figure 7, are significantly better than those of
conventional STPM. Conventional STPM cannot accurately detect anomaly areas of
”Cable”, ”Pill”, and ”Capsule”. In contrast, the proposed method can
accurately detect anomalous regions that STPM cannot. STPM detects the
abnormal area of ”Screw” and ”Transistor,” but there are large anomaly values
in normal areas. The proposed method can accurately detect these anomalies.
These results demonstrate the effectiveness of the proposed method.

TABLE I: Anomaly localization on the MVTec (Pixel level AUROC / PRO AUROC)

| category | PaDim[37] | DRAEM | STPM |  | Ours
---
(Student-Teacher)
| Ours
---
(Discriminative)
Ours
Textures | Carpet | 0.991 / 0.962 | 0.955 / - | 0.988 / 0.958 | 0.988 / 0.952 | 0.985 / 0.952 |  0.992 / 0.968
Grid | 0.973 / 0.946 |  0.997 / - | 0.990 / 0.966 | 0.994 / 0.967 | 0.996 / 0.974 | 0.996 / 0.975
Leather | 0.992 / 0.978 | 0.986 / - | 0.993 / 0.980 | 0.984 / 0.951 |  0.996 / 0.984 |  0.996 / 0.984
Tile | 0.941 / 0.860 |  0.992 / - | 0.974 / 0.921 | 0.970 / 0.910 | 0.987 / 0.958 | 0.988 / 0.962
Wood | 0.949 / 0.911 | 0.964 / - | 0.972 / 0.936 | 0.970 / 0.942 | 0.978 / 0.960 |  0.981 / 0.967
Objects | Bottle | 0.983 / 0.948 | 0.991 / - | 0.988 / 0.951 | 0.990 / 0.965 |  0.993 / 0.962 |  0.993 / 0.973
Cable | 0.967 / 0.888 | 0.943 / - | 0.955 / 0.877 | 0.973 / 0.924 | 0.979 / 0.926 |  0.983 / 0.948
Capsule |  0.985 / 0.935 | 0.943 / - | 0.983 / 0.922 |  0.985 / 0.917 | 0.972 / 0.874 |  0.985 / 0.910
Hazel nut | 0.982 / 0.926 |  0.997 / - | 0.985 / 0.943 | 0.991 / 0.948 | 0.995 / 0.949 | 0.995 / 0.956
Metal nut | 0.972 / 0.856 |  0.995 / - | 0.976 / 0.945 | 0.982 / 0.947 | 0.989 / 0.958 | 0.989 / 0.963
Pill | 0.957 / 0.927 | 0.976 / - | 0.978 / 0.965 | 0.972 / 0.965 |  0.990 / 0.961 | 0.987 / 0.974
Screw | 0.985 / 0.944 | 0.976 / - | 0.983 / 0.930 |  0.993 / 0.958 | 0.989 / 0.935 |  0.993 / 0.955
Toothbrush | 0.988 / 0.931 | 0.981 / - | 0.989 / 0.922 | 0.989 / 0.905 |  0.994 / 0.925 | 0.993 / 0.927
Transistor |  0.975 / 0.845 | 0.909 / - | 0.825 / 0.695 | 0.898 / 0.817 | 0.872 / 0.803 | 0.907 / 0.830
Zipper | 0.985 / 0.959 | 0.988 / - | 0.985 / 0.952 | 0.985 / 0.955 | 0.988 / 0.957 |  0.992 / 0.972
| Mean | 0.975 / 0.921 | 0.973 / - | 0.970 / 0.921 | 0.977 / 0.935 | 0.980 / 0.939 |  0.985 / 0.951

TABLE II: Result of image level AUC on the MVTec

| category | PaDim[37] | DRAEM | STPM |  | Ours
---
(Student-Teacher)
| Ours
---
(Discriminative)
Ours
Textures | Carpet | 0.998 | 0.970 | - | 0.981 | 0.964 | 0.987
Grid | 0.967 | 0.999 | - | 0.984 | 1.000 | 1.000
Leather | 1.000 | 1.000 | - | 0.998 | 1.000 | 1.000
Tile | 0.981 | 0.996 | - | 0.953 | 0.986 | 0.999
Wood | 0.992 | 0.991 | - | 0.992 | 0.966 | 0.993
Objects | Bottle | 0.999 | 0.992 | - | 1.000 | 1.000 | 1.000
Cable | 0.927 | 0.918 | - | 0.967 | 0.989 | 0.996
Capsule | 0.913 | 0.985 | - | 0.873 | 0.945 | 0.930
Hazel nut | 0.920 | 1.000 | - | 1.000 | 0.987 | 0.998
Metal nut | 0.987 | 0.987 | - | 1.000 | 0.999 | 1.000
Pill | 0.933 | 0.989 | - | 0.967 | 0.980 | 0.981
Screw | 0.858 | 0.939 | - | 0.948 | 0.941 | 0.968
Toothbrush | 0.961 | 1.000 | - | 0.900 | 0.987 | 0.979
Transistor | 0.974 | 0.931 | - | 0.975 | 0.978 | 0.983
Zipper | 0.903 | 1.000 | - | 0.898 | 0.998 | 0.993
| Mean | 0.955 | 0.980 | 0.955 | 0.962 | 0.981 | 0.987

TABLE III: Accuracy while changing teacher network for reconstruction student

| ResNet18 | ResNet50(Ours)
---|---|---
Image Level AUC | 0.982 | 0.987
Pixel Level AUC | 0.979 | 0.985
PRO | 0.937 | 0.951

TABLE IV: Accuracy with/without attention mechanism

| Without Attention | All Students–Teachers | Ours
---|---|---|---
Image Level AUC | 0.985 | 0.985 | 0.987
Pixel Level AUC | 0.983 | 0.983 | 0.985
PRO | 0.944 | 0.946 | 0.951

###  IV-C Ablation study

In the proposed method, student network1 and student network2 have different
teacher networks, as shown in Figure 3 and 6, respectively. An attention
mechanism was also used from the teacher network to the student network.
First, the influence of changes in teacher networks is discussed. In this
study, ResNet50 was used as a teacher for student reconstruction. The proposed
method was compared with the method in which ResNet18 was used as a teacher.
Note that all the methods use an attention mechanism from the teacher-to-
student network. Table III lists the accuracy associated with the teacher
network changes. Different teacher networks lead to a perspective different
from that of the conventional method, improving anomaly detection accuracy.

Next, the influence of the attention mechanism is discussed. In this study, an
attention mechanism was used only between Student2 and Teacher2. The proposed
method is compared with the method without the attention mechanism and the
method using attention for both Student1-Teacher1 and Student2-Teacher2. Table
IV shows that it is better to use the attention mechanism only between
Student2 and Teacher2. We believe that the attention mechanism allows us to
learn better. Attention between Student1 and Teacher1 did not improve the
accuracy. Because Student1 and Teacher1 have the same structure, knowledge
distillation is easy. Therefore, if the attention mechanism is used for
student1 and teacher1, the representational ability of the student network is
high, possibly leading to a small difference in features between students and
teachers. The attention mechanism was appropriate for use in Student2 and
Teacher2, providing the expected improvement.

##  V Conclusions

In this study, the STPM anomaly detection method was improved using a new
student network for reconstruction. The teacher network was modified for a new
student network and the attention mechanism from teacher to student was used
to ensure successful learning. To further improve the accuracy of anomaly
detection, discriminative networks were used to reconsider the anomaly maps.
By using pseudo-anomalies in training a discriminative network, a more
accurate anomaly map was generated for normal regions. The proposed method,
using both STPM with two pairs of student-teacher networks and a
discriminative network, achieved high anomaly detection accuracy compared with
conventional methods.

## References

  * [1] Bergmann, P., Fauser, M., Sattlegger, D., Steger, C., Uninformed students: Student-teacher anomaly detection with discriminative latent embeddings. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 4183-4192, 2020.
  * [2] Bergmann, P., Fauser, M., Sattlegger, D., Steger, C., Mvtec-ad real real-world dataset for unsupervised anomaly detection. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 9592-9600, 2019.
  * [3] Bergmann, P., Fauser, M., Sattlegger, D., Steger, C., Uninformed students: Student-teacher anomaly detection with discriminative latent embeddings. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 4183-4192, 2020.
  * [4] Napoletano, P., Piccoli, F., Schettini, R., Anomaly detection in nano brous materials by cnn-based self-similarity. Sensors18(1), 209, 2018.
  * [5] Abati, D., Porrello, A., Calderara, S., Cucchiara, R., Latent space autoregression for novelty detection. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 481-490, 2019.
  * [6] Nguyen, P., Liu, T., Prasad, G., Han, B., Weakly supervised action localization by sparse temporal pooling network. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 6752-6761, 2018.
  * [7] Roitberg, A., Al-Halah, Z., Stiefelhagen, R., Informed democracy: voting-based novelty detection for action recognition. arXiv preprint arXiv:1810.12819, 2018.
  * [8] Baur, C., Wiestler, B., Albarqouni, S., Navab, N., Deep autoencoding models for unsupervised anomaly segmentation in brain MR images. In: International MICCAI Brain lesion Workshop, pp. 161-169, 2018.
  * [9] Li, L., Xu, M., Wang, X., Jiang, L., Liu, H., Attention based glaucoma detection: a large-scale database and CNN model. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 10571-10580, 2019.
  * [10] Schlegl, T., Seeböck, P., Waldstein, S.M., Langs, G., Schmidt-Erfurth, U., f-anoGAN: Fast unsupervised anomaly detection with generative adversarial networks. Medical Image Analysis, 54, 30-44, 2019.
  * [11] Schlegl, T., Seeböck, P., Waldstein, S.M., Schmidt-Erfurth, U., Langs, G., Unsupervised anomaly detection with generative adversarial networks to guide marker discovery. In: International Conference on Information Processing in Medical Imaging, pp. 146-157, 2017.
  * [12] Vasilev, A., Golkov, V., Meissner, M., Lipp, I., Sgarlata, E., Tomassini, V., Jones,D.K., Cremers, D., Q-space novelty detection with variational autoencoders. In: Computational Diffusion MRI, pp. 113-124, 2020.
  * [13] Bergmann, P., Löwe, S., Fauser, M., Sattlegger, D., Steger, C., Improving unsupervised defect segmentation by applying structural similarity to autoencoders. arXiv preprint arXiv:1807.02011, 2018.
  * [14] Wang, G., Han, S., Ding, E., Huang, D., Student-teacher feature pyramid matching for unsupervised anomaly detection. arXiv preprint arXiv:2103.04257, 2021.
  * [15] Nalisnick, E., Matsukawa, A., Teh, Y.W., Gorur, D., Lakshminarayanan, B., Do deep generative models know what they don’t know? arXiv Preprint arXiv:1810.09136, 2018.
  * [16] An, J., Cho, S., Variational autoencoder based anomaly detection using reconstruction probability. Special Lecture on IE 2(1), 1-18, 2015.
  * [17] Zhou, C., Paffenroth, R.C., Anomaly detection with robust deep autoencoders. In: Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 665-674, 2017.
  * [18] Zenati, H., Romain, M., Foo, C.S., Lecouat, B., Chandrasekhar, V., Adversarially learned anomaly detection. In: 2018 IEEE International Conference on Data Mining, pp. 727-736, 2018.
  * [19] Deecke, L., Vandermeulen, R., Ruff, L., Mandt, S., Kloft, M., Image anomaly detection with generative adversarial networks. In: Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pp. 3-17, 2018.
  * [20] Zong, B., Song, Q., Min, M.R., Cheng, W., Lumezanu, C., Cho, D., Chen, H., Deep autoencoding gaussian mixture model for unsupervised anomaly detection. In: International Conference on Learning Representations, 2018.
  * [21] Golan, I., El-Yaniv, R., Deep anomaly detection using geometric transformations. Advances in Neural Information Processing Systems, 31, 2018.
  * [22] Zavrtanik, V., Kristan, M., and Skočaj, D., Reconstruction by inpainting for visual anomaly detection. Pattern Recognition, 112, 107706, 2021.
  * [23] Akcay, S., Atapour-Abarghouei, A., Breckon, T.P., Ganomaly: Semi-supervised anomaly detection via adversarial training. In: Asian Conference on Computer Vision, pp. 622-637, 2018.
  * [24] Akcay, S., Atapour-Abarghouei, A., Breckon, T.P., Skip-ganomaly: Skip connected and adversarially trained encoder-decoder anomaly detection. In: 2019 International Joint Conference on Neural Networks, pp. 1-8, 2019.
  * [25] Tang, T.W., Kuo, W.H., Lan, J.H., Ding, C.F., Hsu, H., Young, H.T., Anomaly detection neural network with dual auto-encoders GAN and its industrial inspection applications. Sensors 20(12), 3336, 2020.
  * [26] Cohen, N., Hoshen, Y., Sub-image anomaly detection with deep pyramid correspondences. arXiv Preprint arXiv:2005.02357, 2020.
  * [27] Wang, F., Jiang, M., Qian C., Yang, S., Li, C., Zhang, H., Wang, X., Tang, X., Residual attention network for image classification. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 3156-3164, 2017.
  * [28] Hu, J., Shen, L., Sun, G., Squeeze-and-excitation networks. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 7132-7141, 2018.
  * [29] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, L., Polosukhin, I., Attention is all you need. In: Advances in Neural Information Processing Systems, pp. 5998-6008, 2017.
  * [30] Fu, J., Liu, J., Tian, H., Li, Y., Bao, Y., Fang, Z., Lu, H., Dual attention network for scene segmentation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 3146-3154, 2019.
  * [31] Ramachandran, P., Parmar, N., Vaswani, A., Bello, I., Levskaya, A., Shlens, J., Stand-alone self-attention in vision models. arXiv preprint arXiv:1906.05909, 2019.
  * [32] Zhang, H., Goodfellow, I., Metaxas, D., Odena, A., Self-attention generative adversarial networks. In: International Conference on Machine Learning, pp.7354-7363, 2019.
  * [33] Fukui, H., Hirakawa, T., Yamashita, T., Fujiyoshi, H., Attention branch network : Learning of attention mechanism for visual explanation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 10705-10714, 2019.
  * [34] Wang, L., Yoon, K.J., Knowledge distillation and student-teacher learning for visual intelligence: A review and new outlooks. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021.
  * [35] Zavrtanik, V., Kristan, M., and Skočaj, D., DRAEM-A discriminatively trained reconstruction embedding for surface anomaly detection. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 8330-8339, 2021.
  * [36] Perlin, K., An image synthesizer. ACM Siggraph Computer Graphics, 19(3), pp. 287-296, 1985.
  * [37] Defard, T., Setkov, A., Loesch, A., and Audigier, R., Padim: a patch distribution modeling framework for anomaly detection and localization. In International Conference on Pattern Recognition, pp. 475-489, 2021.


Conceptually, whenever an event has multiple possible dependency sets, this
composite functor splits such events into multiple distinct events, each with
one dependency set.

In fact the Bruns-Lakser completion of $\mathsf{rdp}(E)$ has a very
interesting interpretation in terms of Merkle trees. Merkle trees were first
introduced as a mechanism for digital signatures [Mer87] but now are used
widely from storage systems to source control to the hash structure of
packages in the Nix package management system. The main idea of the last of
these is that it stores “software components in isolation from each other in a
central component store, under path names that contain cryptographic hashes of
all inputs involved in building the component” [Dol06]. This establishes,
effectively, a tree, where each node can be thought of depending on all nodes
below it, and is labeled by its own data combined with (the recursive hash of)
the subtree of all nodes below it. The act of distinguishing multiple copies
of “the same” event by the record of the choices of their possible inputs
corresponds to the construction of a component store in Nix. There is thus a
sense in which action of the Bruns-Lakser completion on DSCs yields “the
Merkle function”.

##  7\. Version Parametrization

In keeping with our intended notion of DSCs as presenting collections of
packages, this section attempts to give some meaningful semantics to the
informal notion of a “version”. Typically, a version is just some monotone
increasing identifier attached to a computer program or software package with
no special meaning. However, there are often things such as package version
policies, or semantic versioning, which attempt to give some “subtype”-like
meaning to versioning. In particular, the idea is that to upgrade from a lower
version to a higher version should make things in a sense better, rather than
worse – i.e., provide new features, fix bugs, etc. In such version policies,
these universally better, non-breaking updates to a package are known as minor
version releases – as opposed to major version releases, which may change the
functionality arbitrarily, including perhaps scrapping all the existing
features for new ones intended to be used for a similar purpose.

Here we introduce a notion of a version relation that attempts to capture the
idea of such a minor version update, i.e. when one event is universally more
useful than another event. That is to say, when some event is enabled in at
least the situations where the lower version is, and further, enables at least
all the same events enabled by the lower version.

######  Definition 7.1.

Let $(E,\text{dep})$ be a DSC. We introduce a relation on $E$, written
$e\blacktriangleleft e^{\prime}$ if

  1. (1)

$\text{dep}(e)\subseteq\text{dep}(e^{\prime})$, and

  2. (2)

if $D^{x}\in\text{dep}(x)$ is a depset such that $e\in D^{x}$, then
$(D^{x}\setminus e)\cup\\{e^{\prime}\\}\in\text{dep}(x)$ is also a depset of
$x$.

We call this the higher version relation. We say that $e^{\prime}$ is a higher
version of $e$.

######  Lemma 7.2.

The higher version relation $\blacktriangleleft$ is reflexive and transitive.

###### Proof.

It is clear that $\blacktriangleleft$ is reflexive. Now suppose that
$e\blacktriangleleft e^{\prime}$ and $e^{\prime}\blacktriangleleft
e^{\prime\prime}$. Then if $D^{x}\in\text{dep}(x)$ for some $x\in E$ with
$e\in X$, then $D^{x}_{0}=(D^{x}\setminus
e)\cup\\{e^{\prime}\\}\in\text{dep}(x)$. But since
$e^{\prime}\blacktriangleleft e^{\prime\prime}$, and $e^{\prime}\in
D^{x}_{0}$, then
$(D^{x}_{0}\setminus\\{e^{\prime}\\})\cup\\{e^{\prime\prime}\\}=(D^{x}\setminus
e)\cup\\{e^{\prime\prime}\\}\in\text{dep}(x)$. Thus $e\blacktriangleleft
e^{\prime\prime}$. ∎

Given $e\in E$, let $\text{Vers}(e)$ denote the set of $e^{\prime}\in E$ such
that $e\blacktriangleleft e^{\prime}$. We call this the set of higher versions
of $e$.

Now if $(E,\text{dep})$ is a DSC, then consider the following function,
$V:P(E)\to P(E)$ defined by

$V(X)=\bigcup_{x\in X}\text{Vers}(x)$

######  Lemma 7.3.

If $X\subseteq E$ is a reachable subset, then $V(X)$ is a reachable subset.

###### Proof.

Suppose that $e^{\prime}\in V(X)$. Then there exists some $e\in X$ such that
$e\blacktriangleleft e^{\prime}$. Then
$\text{dep}(e)\subseteq\text{dep}(e^{\prime})$. Since $X$ is complete, there
exists some $D^{e}\subseteq X$. But $D^{e}\in\text{dep}(e^{\prime})$ and
$D^{e}\subseteq X\subseteq V(X)$, thus $V(X)$ is complete and therefore
reachable. ∎

######  Definition 7.4.

Let $P$ be a poset and $f:P\to P$ a function. We say that $f$ is a closure
operator if:

  1. (1)

$f$ is order preserving,

  2. (2)

if $X\in P$, then $X\leq f(X)$, and

  3. (3)

if $X\in P$, then $f(f(X))=f(X)$.

Note that this is precisely the same thing as an idempotent monad on $P$
thought of as a category.

######  Corollary 7.5.

Given a DSC $(E,\text{dep})$, $V$ descends to a function
$V:\mathsf{rdp}(E)\to\mathsf{rdp}(E)$, and is a closure operator.

Now we wish to describe how $V$ descends to a closure operator on the Bruns-
Lakser completion.

######  Lemma 7.6.

Let $L$ denote a finite lattice, and suppose that
$S\in\mathcal{O}\mathcal{J}(L)$. Then

$S=\bigcup_{x\in S}\varphi_{L}(x).$

###### Proof.

Suppose that $x\in S$. Then $x$ is join irreducible, so $x\in\varphi_{L}(x)$,
thus $x\in\bigcup_{x\in S}\varphi_{L}(x)$ and so $S\subseteq\bigcup_{x\in
S}\varphi_{L}(x)$.

Now suppose that $y\in\bigcup_{x\in S}\varphi_{L}(x)$. Then
$y\in\varphi_{L}(x)$ for some $x\in S$. Thus $y\leq x$ and $y$ is join
irreducible. But $S$ is a down-closed subset of $\mathcal{J}(L)$, so $y\in S$,
and therefore $\bigcup_{x\in S}\varphi_{L}(x)\subseteq S$. ∎

######  Corollary 7.7.

Let $L$ denote a finite lattice and suppose $S\in\mathsf{BL}(L)$, then

$S=\bigcup_{x\in S}\mathbb{b}(x).$

###### Proof.

By Proposition 6.14. ∎

If $V:L\to L$ is a closure operator on a finite lattice, we want to show that
it induces a closure operator $\mathbb{V}:\mathsf{BL}(L)\to\mathsf{BL}(L)$. If
$S\in\mathsf{BL}(L)$, then by Corollary 7.7, we can write $S$ as
$S=\bigcup_{x\in S}\mathbb{b}(x)$. Thus define

$\mathbb{V}(S)=\bigcup_{x\in S}\mathbb{b}(V(x)).$

######  Proposition 7.8.

The function $\mathbb{V}:\mathsf{BL}(L)\to\mathsf{BL}(L)$ as defined above is
a closure operator.

###### Proof.

Suppose that $S\subseteq T$. Then clearly
$\mathbb{V}(S)\subseteq\mathbb{V}(T)$. Thus $\mathbb{V}$ is order preserving.

We wish to show that $S\subseteq\mathbb{V}(S)$. Suppose that $x\in S$. Then
$x\leq V(x)$, since $V$ is a closure operator on $L$. Thus
$\mathbb{b}(x)\subseteq\mathbb{b}(V(x))$. Therefore $S\subseteq\bigcup_{x\in
S}\mathbb{b}(V(x))=\mathbb{V}(S)$.

Now we wish to show that $\mathbb{V}(\mathbb{V}(S))\subseteq\mathbb{V}(S)$.
Suppose that $z\in\mathbb{V}(\mathbb{V}(S))$. Then $z\in\mathbb{b}(V(y))$ for
some $y\in\mathbb{V}(S)$. But then $y\in\mathbb{b}(V(x))$ for some $x\in S$.
This implies $y\leq V(x)$, so $V(y)\leq V(V(x))\leq V(x)$, since $V$ is a
closure operator. Thus $\mathbb{b}(V(y))\subseteq\mathbb{b}(V(x))$, and
$z\in\mathbb{b}(V(x))\subseteq\mathbb{V}(S)$, therefore $\mathbb{V}$ defines a
closure operator. ∎

Thanks to Proposition 7.8, the closure operator
$\mathbb{V}:\mathsf{BL}(\mathsf{rdp}(E))\to\mathsf{BL}(\mathsf{rdp}(E))$ is
now defined for any DSC $(E,\text{dep})$. We call this the version monad on
$\mathsf{BL}(\mathsf{rdp}(E))$.

##  8\. Conclusion

From a mathematical standpoint, this paper has made the following
contributions. Theorem 3.18 shows that antimatroids can be used to model
package management systems. By introducing DSCs we have provided a connection
between the theory of general event structures and that of antimatroids. We
have also proposed a natural definition of morphisms of antimatroids, and
explored properties of the resulting category of antimatroids. Furthermore, we
have given and shown applications of a simple characterization of the finite
Bruns-Lakser completion.

From an applied standpoint, we have put this mathematics to work in providing
a formal characterization of the dependency data of a package repository, as
well as a formal characterization of “minor” versions of packages, and of
package versioning policies. We have also given a mathematical account of the
logical and order-theoretic interpretation of Merkle trees, which sheds light
on the reasons for their widespread utility.

It is worth noting that the construction and characterization of Bruns-Lakser
given in Proposition 6.14 in fact extends beyond finite lattices to all finite
posets, as explored in [BP20].

There remain some interesting outstanding questions regarding morphisms of
DSCs. It was established in Corollary 2.42 that the image of morphisms of DSCs
under $\mathsf{rdp}$ is a map of join-semilattices. It is an open question if
this is a bi-implication. Similarly, it is open if Lemma 2.46 is a bi-
implication, and if not, what an appropriate collection of morphisms would be
to generate an equivalence of categories between some subcategory of
$\mathsf{DSC}$ and finite lattices. Finally, it is an open question to
characterize a larger class of maps of DSCs that makes $\mathsf{BL}$
functorial.

There are a number of possible directions for future work. One is further
development of a category whose internal logic can encode package dependency
problems. Another is modeling updates to package repositories by means of
morphisms in the category of DSCs. Still another is extending DSCs with a
notion of conflict, so that they would more fully correspond to general event
structures. Beyond that is development of insights into efficient dependency
solvers – a topic of great practical interest. Finally, there is much we hope
that can be extracted from the surprising connection between dependency
systems and antimatroids. There are many high-powered tools available for the
study of the latter, as well as their cousins, matroids, and the more general
class of greedoids to which both belong, not least characteristic polynomials.
It would be of interest to see what insights these tools can bring when
transported to the world of dependency structures.

## References

  * [Baz21] Gershom Bazerman  “The Semantics of Package Management via Event Structures”, 2021  arXiv:2107.01542
  * [BF90] E Andrew Boyd and Ulrich Faigle  “An algorithmic characterization of antimatroids”  In _Discrete Applied Mathematics_ 28.3 Elsevier, 1990, pp. 197–205
  * [Bir40] Garrett Birkhoff  “Lattice theory”  American Mathematical Soc., 1940
  * [BL70] Günter Bruns and Harry Lakser  “Injective hulls of semilattices”  In _Canadian Mathematical Bulletin_ 13.1 Cambridge University Press, 1970, pp. 115–118
  * [BP20] Gershom Bazerman and Raymond Puzio  “The Topological and Logical Structure of Concurrency and Dependency via Distributive Lattices”, 2020  arXiv:2004.05688
  * [Czé14] Gábor Czédli  “Coordinatization of finite join-distributive lattices”  In _Algebra universalis_ 71.4 Springer ScienceBusiness Media LLC, 2014, pp. 385–404
  * [Ded00] Richard Dedekind  “Über die von drei Moduln erzeugte Dualgruppe”  In _Mathematische Annalen_ 53.3 Springer, 1900, pp. 371–403
  * [Dil40] R.. Dilworth  “Lattices with Unique Irreducible Decompositions”  In _Annals of Mathematics_ 41.4 Annals of Mathematics, 1940, pp. 771–777
  * [Dol06] Eelco Dolstra  “The purely functional software deployment model”  Utrecht University, 2006
  * [Ede80] Paul H Edelman  “Meet-distributive lattices and the anti-exchange closure”  In _Algebra Universalis_ 10.1 Springer, 1980, pp. 290–299
  * [GG13] Mai Gehrke and Sam Van Gool  “Distributive envelopes and topological duality for lattices via canonical extensions”, 2013  arXiv:1309.3113
  * [GY94] Paul Glasserman and David D Yao  “Monotone structure in discrete-event systems”  John Wiley & Sons, Inc., 1994
  * [HP18] Chris Heunen and Vaia Patta  “The category of matroids”  In _Applied Categorical Structures_ 26.2 Springer, 2018, pp. 205–237
  * [Joh82] Peter T Johnstone  “Stone spaces”  Cambridge university press, 1982
  * [KL13] Yulia Kempner and Vadim E Levit  “Geometry of poset antimatroids”  In _Electronic Notes in Discrete Mathematics_ 40 Elsevier, 2013, pp. 169–173
  * [KL85] Bernhard Korte and László Lovász  “Shelling structures, convexity and a happy end”  Inst. für Ökonometrie u. Operations-Research, 1985
  * [Mer87] Ralph C Merkle  “A digital signature based on a conventional encryption function”  In _Conference on the theory and application of cryptographic techniques_ , 1987, pp. 369–378  Springer
  * [Mon+18] Ugo Montanari et al.  “Event Structures for Petri nets with Persistence”  In _Logical Methods in Computer Science_ 14 Episciences. org, 2018
  * [Mon85] Bernard Monjardet  “A use for frequently rediscovering a concept”  In _Order_ 1.4 Springer, 1985, pp. 415–417
  * [MR16] Nazarré Merchant and Jason Riggle  “OT grammars, beyond partial orders: ERC sets and antimatroids”  In _Natural Language & Linguistic Theory_ 34.1 Springer, 2016, pp. 241–269
  * [NPW79] Mogens Nielsen, Gordon Plotkin and Glynn Winskel  “Petri nets, event structures and domains”  In _Semantics of Concurrent Computation_ Springer, 1979, pp. 266–284
  * [Par03] Aarati Parmar  “Some mathematical structures underlying efficient planning”  In _AAAI spring symposium on logical formalization of commonsense reasoning_ , 2003
  * [Stu05] Isar Stubbe  “The canonical topology on a meet-semilattice”  In _International Journal of Theoretical Physics_ 44.12 Springer, 2005, pp. 2283–2293
  * [XS15] Zhen-Yu Xiu and Fu-Gui Shi  “Categories of-Fuzzy Greedoids”  In _Journal of Function Spaces_ 2015 Hindawi, 2015


# Plotting in a Formally Verified Way

Guillaume Melquiond This work was funded by the NuSCAP project
(ANR-20-CE48-0014) of the French national research agency (ANR).Université
Paris-Saclay, CNRS, ENS Paris-Saclay, Inria, Laboratoire Méthodes Formelles,
Gif-sur-Yvette, 91190, France<EMAIL_ADDRESS>
###### Abstract

An invaluable feature of computer algebra systems is their ability to plot the
graph of functions. Unfortunately, when one is trying to design a library of
mathematical functions, this feature often falls short, producing incorrect
and potentially misleading plots, due to accuracy issues inherent to this use
case. This paper investigates what it means for a plot to be correct and how
to formally verify this property. The Coq proof assistant is then turned into
a tool for plotting function graphs using reliable polynomial approximations.
This feature is provided as part of the CoqInterval library.

##  1 Introduction

An invaluable feature of computer algebra systems (Maple, Mathematica, etc) is
their ability to plot the graph of mathematical functions. Indeed, as the
adage goes, a picture is worth a thousand words. When encountering an unknown
function, the first reflex of the user is to plot it, so as to grasp its
features. But how much can the user trust that the plot is an accurate
depiction of the graph of the function?

Let us consider the case of a user who wants to implement a floating-point
library of mathematical functions. The implementation of such a function,
_e.g._ , $\exp$, usually involves some polynomial $p$, since processors are
efficient at addition and multiplication. The distance $|p(x)-\exp x|$
characterizes the quality of the approximation and thus of the implementation.
So, one might be tempted to plot $p(x)-\exp x$ and look at its extrema. Let us
assume that $p$ is the minimax approximation of degree 6 between $-2^{-5}$ and
$2^{-5}$ with _binary64_ floating-point coefficients. Figure 1(a) shows what
the signed distance looks like, when plotted with Gnuplot. The result
certainly looks questionable, and other computer algebra systems would hardly
do better.

The main issue is that these systems plot the function graph using plain
_binary64_ arithmetic, which is way too inaccurate for our use case. The
Sollya tool was especially designed to solve this kind of issue [2]. By using
a 165-bit arithmetic, it is able to produce Figure 1(c), which is
representative of what the distance between a function and its minimax
polynomial usually looks like (_cf._ La Vallée-Poussin’s theorem). The user
can now look at the plotted function graph and see that the distance is
bounded by $10^{-16}$, which might be sufficient, depending on the purpose of
the mathematical library.

(a) Gnuplot

(b) Sollya

(c) Sollya and this work

(d) This work

Figure 1: On the left, distance between $\exp$ and its minimax. On the right,
$x\to\sin(x+\exp x)$.

But what if the function graph needs more than 165 bits of precision to be
plotted correctly? Sollya actually performs its computations using interval
arithmetic. Instead of computing a single floating-point number that
approximates the real value of the function, it computes two floating-point
bounds that enclose this real value. So, by looking at the distance between
these bounds, it is able to detect when its internal precision is not
sufficient for a correct plot. In that case, the user can increase the
precision.

So, did Sollya actually solve the issue of plotting? Not quite. As other
computer algebra systems, it heuristically samples the plotted function. So,
any feature of the function that occurs between two sampled abscissas will not
appear on the plot. Consider the function $f(x)=\sin(x+\exp x)$. This time,
there is no precision issue; even _binary32_ floating-point arithmetic could
be used. But the function oscillates so quickly that the plot is not faithful,
as shown on Figure 1(b). For $x\geq 4$, the function does not even seem to
reach $-1$ and $1$ anymore, although it is the sine function.

Let us define a _correct plot_ as follows: If a pixel is blank, then there
exists no value of $x$ such that $(x,f(x))$ falls into this pixel. Conversely,
a _complete plot_ is defined as follows: If a pixel is filled, then there
exists a value of $x$ such that $(x,f(x))$ falls into it. This article
explains how one can draw correct plots, as shown on Figures 1(c) and 1(d). To
increase the confidence in the plots, they are performed using the Coq system.
This is not the first work to explore the topic of correct (and complete)
plots using a proof assistant [6], but this works focuses on producing plots
in a matter of seconds rather than hours.111Figure 1(d) needs more than
600,000 samples when using a method based on the modulus of uniform continuity
[6].

Section 2 explains how to formally define what a correct plot is. Section 3
then shows how to compute it inside the logic of Coq, using the CoqInterval
library [5]. It also explains how to get close to a complete plot. Coq gives
the greatest confidence in the correctness of the plots, but it hardly strikes
as a user-friendly system when it comes to computer algebra. So, Section 4
focuses on interface concerns.

##  2 Formal Specification

The very first step is to state what it means for a plot to be correct. We
have some function $f$ from real numbers to real numbers. Note that this is a
function defined in a purely mathematical sense; no floating-point numbers are
involved here. We also have some bounds $x_{\\_}1$ and $x_{\\_}2$ between
which we will plot the function graph. Again, these are real numbers. We will
not use them directly though, as it makes for much more readable (and thus
trustable) definitions to use two real numbers $\mathit{ox}$ and $\mathit{dx}$
such that $x_{\\_}1=\mathit{ox}$ and $x_{\\_}2=\mathit{ox}+\mathit{dx}\cdot
w$. In other words, if the output device has an horizontal resolution of $w$
pixels, then $\mathit{dx}$ will be the width of a single pixel, while
$\mathit{ox}$ will point at the left border of the leftmost pixel.

A plot will then be defined as a list $\ell$ of intervals. The $i$-th interval
$\ell_{\\_}i$ encloses all the possible values of the function for the $i$-th
pixel:

$\forall x,~{}\mathit{ox}+\mathit{dx}\cdot i\leq
x\leq\mathit{ox}+\mathit{dx}\cdot(i+1)\Rightarrow f(x)\in\ell_{\\_}i.$

The translation to Coq is straightforward (with I a module containing
interval-related definitions):

⬇

Definition plot1 (f:R->R) (ox dx:R) (l:list I.type) :=

forall (i:nat) (x:R), ox + dx * i <= x <= ox + dx * (i + 1) ->

I.contains (nth i l I.nai) (f x).

The interval I.nai, which contains any real number, is the default value used
when the list $\ell$ is exhausted. Thus, the predicate plot1 is still valid
when $i$ exceeds the length of $\ell$. In particular, if this length is
smaller than $w$, the rightmost part of the plot will be entirely comprised of
filled pixels.

This could be the end of it, but since I.type is an internal datatype, it
might be difficult to turn the list into a bitmap or to serialize it to
another tool. So, in order to use a more universal datatype, namely integers
of type Z, a second predicate is defined. To do so, we need two more real
numbers: $\mathit{oy}$ and $\mathit{dy}$. They play a role similar to
$\mathit{ox}$ and $\mathit{dx}$, but along the vertical axis:

⬇

Definition plot2 (f:R->R) (ox dx oy dy:R) (h:Z) (l:list (Z*Z)) :=

forall (i:nat) (x:R), ox + dx * i <= x <= ox + dx * (i + 1) ->

oy <= f x <= oy + dy * h ->

let r := nth i l (0, h) in

oy + dy * (fst r) <= f x <= oy + dy * (snd r).

The seemingly counterproductive hypothesis $\mathit{oy}\leq
f(x)\leq\mathit{oy}+\mathit{dy}\cdot h$ is there in case the user wants to
focus on some detail of the function graph and is not interested in extreme
values that happen near it. As translating these extreme values to integers is
useless and possibly dangerous (due to overflow), only the pixels with
ordinates between $0$ and $h$ are kept, with $h$ the vertical resolution of
the output device.

If the user did not explicitly provide $y_{\\_}1=\mathit{oy}$ and
$y_{\\_}2=\mathit{oy}+\mathit{dy}\cdot h$, they can be found by computing the
union of all the intervals of a list satisfying plot1. In that case, the
hypothesis $\mathit{oy}\leq f(x)\leq\mathit{oy}+\mathit{dy}\cdot h$ is
trivially satisfied. This explains why the proposed mechanism involves two
predicates plot1 and plot2, instead of having just plot2. The following is an
instance obtained for the function $x\mapsto x^{2}$ between $0$ and $1$ for
$w=10$ and $h=100$. Notice how $\mathit{oy}\simeq-0.0003\leq 0$ and
$\mathit{oy}+\mathit{dy}\cdot h\simeq 1.015\geq 1$.

⬇

plot2 (fun x => x^2) 0 (820/8192) (-5/16384) (665/65536) 100

((0, 2) :: (0, 5) :: (3, 9) :: (8, 16) :: (15, 25) :: (24, 36)

:: (35, 49) :: (48, 64) :: (62, 81) :: (79, 100) :: nil)

##  3 Plotting a function graph

Converting a list that satisfies plot1 into a list that satisfies plot2 is a
bit technical, but there is no difficulty in formally verifying the algorithm
using Coq. Thus, let us focus on getting a list of intervals that satisfies
plot1.

The first step is to turn the function $f$ into a straight-line program that
can be manipulated by functions written in Gallina, the language of Coq, a
dependently-typed lambda-calculus with inductive datatypes. To do so, we just
reuse the machinery of the interval tactic, that is, an Ltac oracle reifies
the function and Coq then formally checks that the reified function matches
the original one [5]. For example, if the original function is $x\mapsto\cos
x+3$, the reified function will look like (Binary Add (Unary Cos (Var 0))
(Const 3)) where Binary, Cos, etc are constructors of some inductive
datatypes.

The second step is to compute the list $\ell$ such that
$f(X_{\\_}i)\subseteq\ell_{\\_}i$ for $X_{\\_}i=[\mathit{ox}+\mathit{dx}\cdot
i;\mathit{ox}+\mathit{dx}\cdot(i+1)]$. The natural idea would be to use
interval arithmetic to do so. Indeed, for every operation $\diamond$ over the
real numbers, it provides an operation over intervals (abusively noted
$\diamond$ too) that is compatible with it:

$\forall U,V\in\mathbb{I},~{}\forall u,v\in\mathbb{R},~{}u\in U\land v\in
V\Rightarrow u\diamond v\in U\diamond V.$

So, we would just have to recursively visit the reified expression, applying
the corresponding interval operations along the way. This would give us some
interval $\ell_{\\_}i$ that satisfies plot1. While it might work most of the
times, it is not suitable for the use case described in the introduction.
Indeed, the core defect of interval arithmetic, _i.e._ , loss of correlation,
applies here. A pixel might be too wide for interval arithmetic to produce a
meaningful interval. In other words, $\ell_{\\_}i$ will contain $f(x)$, but it
will also contain many other ordinates. Figure 2(a) shows the result; not only
is it a massive block of pixels, but the bounds are off by several orders of
magnitude: $6\cdot 10^{-5}$ instead of $10^{-16}$.

(a) Naive interval arithmetic.

(b) Completeness check disabled.

Figure 2: Correct but hardly complete variants of Figure 1(c).

There is a second, more practical, reason for not using naive interval
arithmetic. Since we are performing all our computations in the logic of Coq,
evaluating $w=512$ instances of the interval implementation of cosine might
take too long for an interactive use of Coq. The solution is to compute a
rigorous polynomial approximation $(p,\Delta)$ of $f$ over some large interval
$X$, ideally $W=[\mathit{ox};\mathit{ox}+\mathit{dx}\cdot w]$:

$\forall x\in X,~{}p(x)-f(x)\in\Delta.$

The idea of using these polynomial approximations originated from Sollya [2].
They were later formalized in Coq [4]. Eventually, they joined the CoqInterval
library [5].

These polynomial approximations were instrumental when devising the integral
tactic for guaranteed numerical quadrature [3]. Indeed, once $(p,\Delta)$ has
been computed, one can easily enclose the integrals of $p(x)$ and of
$p(x)-f(x)$ over $X$, and thus the integral of $f(x)=p(x)-(p(x)-f(x))$ over
$X$. Since the ability to numerically integrate a function is not that
different from the ability to accurately plot its graph, we follow a similar
approach here. Once $(p,\Delta)$ has been computed, we use it to compute an
enclosure $Y_{\\_}i$ of $f(X_{\\_}i)$. This time, we can use naive interval
arithmetic to evaluate $Y_{\\_}i=p(X_{\\_}i)+\Delta$.

If $p$ has degree $0$, then the result is similar to the one obtained using
naive interval arithmetic. With degrees $1$ and $2$, the plot is still an
indiscriminate block of pixels, though the bounds are less overapproximated.
With degree $3$, losses of correlation at the pixel level are completely
accounted for. So, the plot looks fine, but computing it is way too slow for
interactive use. Indeed, a single polynomial is not sufficient for the whole
interval $W$, since there is no way a degree-3 polynomial could ever
meaningfully approximate a function with 7 roots. So, the interval $W$ has to
be split into many subintervals $X$, on each of which $p$ has to be computed.
The best running time is obtained for degree 6. Then, the higher the degree,
the slower it gets. Indeed, decreasing the number of subintervals no longer
compensates the increasing cost of computing $p$ and evaluating it for every
pixel. The optimal degree highly depends on the plot, so its choice is left to
the user. By default, degree 10 is used, as with tactic integral.

Thanks to the rigorous polynomial approximations, we now have a correct plot
that is formally verified. But as shown on Figure 2(a), a correct plot is not
necessarily a complete one. So, to increase the usability of our approach, we
would like the plot to never be more than a few pixels wide. A first idea
would be to measure the width of $p(X_{\\_}i)$. If it is larger than a few
pixels, then $X$ needs to be subdivided further. Unfortunately, this causes
too many subdivisions when the function varies quickly, which is the case at
the left and right ends of Figure 1(c). A second idea would be to measure the
width of the error interval $\Delta$ to decide whether $X$ is sufficiently
small. Unfortunately, this still does not work. Indeed, the further from the
center of $X$, the worst the loss of correlation becomes when evaluating
$p(X_{\\_}i)$, to the point where it becomes noticeable, as shown on Figure
2(b). So, this time, there are not enough subdivisions.

To strike a balance between these two issues, the code computes an
underestimation of $f$ over $X_{\\_}i$ and compares it to the overestimation
$Y_{\\_}i=p(X_{\\_}i)+\Delta$. If the latter is only a few pixels larger than
the former, then it is deemed good enough. Concretely, the code computes
$Z=p(\mathrm{lower}(X_{\\_}i))+\Delta$. There is a value of $x\in X_{\\_}i$
such that $f(x)$ lies between the lower bound of $Y_{\\_}i$ and the upper
bound of $Z$. So, if the distance between these two bounds is smaller than a
few pixels, the lower bound of $Y_{\\_}i$ is accurate enough. If not, the code
tries again with the upper bound of
$Z^{\prime}=p(\mathrm{upper}(X_{\\_}i))+\Delta$. If the lower bound of
$Y_{\\_}i$ is accurate enough, the code then checks the upper bound by
comparing it to the lower bounds of $Z$ and $Z^{\prime}$.

To finish, there is some kind of a chicken-or-egg problem. If the user has not
provided $\mathit{dy}$, how does the code know the height of a pixel used to
check for pseudo-completeness? To estimate it, the code samples the function
at 50 uniformly spaced points of $W$. This gives an underestimation of
$[y_{\\_}1;y_{\\_}2]$ and thus of $\mathit{dy}$. If the sampled values do not
capture the extreme values of the function, then the predicted value of
$\mathit{dy}$ is too small, which causes the plot to be uselessly accurate and
thus slower to compute.

##  4 Interface

We now have an algorithm (run inside the logic of Coq) that, given some
reification of function $f$ and some values for $\mathit{ox}$, $\mathit{dx}$,
etc, computes a list $\ell$ of pairs of integers. We also have a theorem
formally verified in Coq that states (plot2 $f$ $\mathit{ox}$ … $\ell$). So,
only interface issues remain.

First, let us deal with the plot display. The list $\ell$ is almost a run-
length encoding of the function graph, assuming a column-major order. Indeed,
a pair $(y_{\\_}1,y_{\\_}2)\subseteq[0;h]$ of integers represents a column of
first $y_{\\_}1$ blank pixels, then $y_{\\_}2-y_{\\_}1$ filled pixels, and
finally $h-y_{\\_}2$ blank pixels. Getting Gnuplot to draw the resulting
bitmap is easy. Unfortunately, it is difficult to make sure that Gnuplot maps
one pixel of the bitmap to exactly one pixel of the screen. As a consequence,
some features of the plot might disappear if the drawing area is just one
pixel too small. Conversely, if the user tells Gnuplot to zoom in (or just
enlarges the drawing window), then the plot starts looking blocky.

So, rather than a bitmap, it is visually more satisfying to turn the plot into
two piecewise affine curves that enclose the filled pixels of the bitmap. This
vector encoding allows the user to freely zoom on the plot or resize the
Gnuplot windows. Computing these two curves is actually quite easy. Given two
consecutive elements of the list $\ell_{\\_}{i-1}=(y_{\\_}1,y_{\\_}2)$ and
$\ell_{\\_}i=(y^{\prime}_{\\_}1,y^{\prime}_{\\_}2)$, one just needs to
associate to abscissa $\mathit{ox}+i\cdot\mathit{dx}$ the ordinates
$\mathit{oy}+\min(y_{\\_}1,y^{\prime}_{\\_}1)\cdot\mathit{dy}$ and
$\mathit{oy}+\max(y_{\\_}2,y^{\prime}_{\\_}2)\cdot\mathit{dy}$. The band
between the two curves contains all the filled pixels of the original bitmap,
thus guaranteeing the correctness of the plot, at the expense of being a bit
less narrow than the bitmap one.

As for the user queries, let us take some inspiration from existing computer
algebras system. They often handle plots as first-class citizen, that is, the
user can execute “p := plot(f,x1,x2)” to store a function graph into some
variable p. Then, simply executing p causes the plot to be displayed. (Both
steps can usually be merged into a single one, if the function graph does not
need to be stored for later use.) We can follow a similar approach for Coq.
Unfortunately, Coq requires top-level terms to be preceded by a command,
_e.g._ , Print, Check, About. We cannot reuse an existing command, so we add
yet another one: Plot p. This causes Coq to open a Gnuplot windows using the
data encoded in the plot2 type of p. As for a Coq equivalent to “p :=
plot(f,x1,x2)”, we combine the Definition command with the tactic-in-term
feature of Coq. This provides the following interface: Definition p :=
ltac:(plot f x1 x2). The plot tactic can take two extra arguments to specify
the ordinate range. When absent, the tactic computes the extrema of the
function between the endpoints.

The plot tactic supports the same configuration mechanism as the tactics
interval and integral [5]. For instance, to obtain Figure 1(c), one needs to
increase the precision to 90 bits, as follows:

⬇

plot (fun x => 1+x*... - exp x) (-1/32) (1/32) with (i_prec 90)

A new flag has been added to specify the dimension of the bitmap: i_size w h.
By default, the tactic produces a plot of size $512\times 384$. Other
meaningful flags are i_degree to control the degree of the polynomial
approximations and i_native_compute to tell Coq to first compile the algorithm
to machine code rather than directly interpreting it. This might be useful to
speed up some computationally-intensive plots. Indeed, the architecture of Coq
unfortunately forces the tactic to execute the algorithm twice: once to get
the actual list and a second time to instantiate the correctness theorem.

##  5 Conclusion

This article has presented a mechanism integrated in release 4.2 of the
CoqInterval222https://coqinterval.gitlabpages.inria.fr/ library. It makes it
possible to compute formally correct function graphs and display them directly
from Coq. Despite the computations being performed inside the logic of Coq,
performances are good enough for interactive use. For example, it takes less
than 4 seconds to compute and formally verify the complicated plot of Figure
1(d) with the default settings, and about 1 second with i_native_compute.

While the current interface is a bit unfriendly, we could readily imagine a
new front-end that would exempt the user from typing Definition and Plot, as
well as the ltac:(...) quotation mechanism. The long-term goal is to revisit
the way Coq is used, making it more of a computer algebra system, _e.g._ ,
through interfaces such as CoCalc and Jupyter [7]. Plots would no longer be
opened in separate windows but directly embedded in the document. In the
meantime, by virtue of the plotting algorithm being written in Gallina, it
could easily be extracted to OCaml and distributed as a standalone library.

## References

  * [1]
  * [2] Sylvain Chevillard, Mioara Joldeş & Christoph Lauter (2010): _Sollya: An Environment for the Development of Numerical Codes_.  In Komei Fukuda, Joris van der Hoeven, Michael Joswig & Nobuki Takayama, editors: 3rd International Congress on Mathematical Software (ICMS), Lecture Notes in Computer Science 6327, Kobe, Japan, pp. 28–31, 10.1007/978-3-642-15582-65.
  * [3] Assia Mahboubi, Guillaume Melquiond & Thomas Sibut-Pinote (2019): _Formally Verified Approximations of Definite Integrals_.  Journal of Automated Reasoning 62(2), pp. 281–300, 10.1007/s10817-018-9463-7.
  * [4] Érik Martin-Dorel, Micaela Mayero, Ioana Pasca, Laurence Rideau & Laurent Théry (2013): _Certified, Efficient and Sharp Univariate Taylor Models in Coq_.  In: 15th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC), Timisoara, Romania, pp. 193–200, 10.1109/SYNASC.2013.33.
  * [5] Érik Martin-Dorel & Guillaume Melquiond (2016): _Proving Tight Bounds on Univariate Expressions with Elementary Functions in Coq_.  Journal of Automated Reasoning 57(3), pp. 187–217, 10.1007/s10817-015-9350-4.
  * [6] Russell O’Connor (2008): _A Computer Verified Theory of Compact Sets_.  In: Symbolic Computation in Software Science Austrian-Japanese Workshop (SCSS), RISC-Linz Report Series, pp. 148–162.
  * [7] Fernando Pérez & Brian E. Granger (2007): _IPython: a System for Interactive Scientific Computing_.  Computing in Science and Engineering 9(3), pp. 21–29, 10.1109/MCSE.2007.53.