# A Sui Generis QA Approach using RoBERTa for Adverse Drug Event
Identification

Harshit Jain<EMAIL_ADDRESS>ZS AssociatesBengaluruIndia ,  Nishant
Raj<EMAIL_ADDRESS>ZS AssociatesBengaluruIndia and  Suyash Mishra
<EMAIL_ADDRESS>ZS AssociatesLondonUnited Kingdom

###### Abstract.

Extraction of adverse drug events from biomedical literature and other textual
data is an important component to monitor drug-safety and this has attracted
attention of many researchers in healthcare. Existing works are more pivoted
around entity-relation extraction using bidirectional long short term memory
networks (Bi-LSTM) which does not attain the best feature representations. In
this paper, we introduce a question answering framework that exploits the
robustness, masking and dynamic attention capabilities of RoBERTa by a
technique of domain adaptation and attempt to overcome the aforementioned
limitations. Our model outperforms the prior work by 9.53% F1-Score.

adverse drug event, RoBERTa, question-answering, entity-relation extraction,
healthcare

##  1\. Introduction

Recent advancements in drug development and approval have bolstered our
healthcare ecosystem. However, this rapid paced development is accompanied by
an increase in associated risks. Adverse Drug Events (ADEs) form an integral
component of those risks. An ADE is defined as ”an injury resulting from a
medical intervention related to a drug”(Bates et al., 1995).

These events create an economic burden over the system. A study by
Rocchiccioli et al. showed ”a statistically significant increase in all direct
costs (inpatient,outpatient, therapy) during post-ADE period (+US$1310 for all
ADEs and +US$1983 for preventable ADEs, versus the pre-event period)”
(Rocchiccioli et al., 2007; Wu et al., 2012; Chiatti et al., 2012). Further,
national estimates are indicative of the fact that ADEs contribute at least an
additional US$ 30 billion to US healthcare costs (Jagannatha et al., 2019;
Johnson and Booman, 1996). Thus, it becomes an utmost requirement to identify
these ADEs at an early stage from biomedical literatures, EHR data and other
sources in order to avoid additional costs incurred in patient management
while improving pharmacovigilence practices at the same time.

Several Adverse Drug Events have been reported to U.S. Food & Drug
Administration (FDA) through Federal Adverse Event Reporting System (FAERS).
These reports submissions are voluntary and hence they may suffer from cases
of massive under-reporting (Edlavitch, 1988). Hence, researchers have started
moving towards more automated approaches in machine learning. There has been a
gradual shift towards using natural language processing (NLP) based methods.
Early attempts have incorporated the use of resources like NLM’s MetaMap,
Unified Medical Language System (UMLS) etc. (Zelenko et al., 2003) to extract
drugs for ADE identification tasks. However, a major limitation of these
approaches is that they are not able to capture the causal relationships
between drug and ADE properly.

Another popular approach in the community has been to treat this problem as
entity recognition and relation identification task. This approach has further
been tackled using either a 1\. pipeline method where entity recognition tasks
are done first followed by relation identification or 2\. a joint training
method where training weights are shared between these two sub-tasks so that
errors don’t accumulate (Miwa and Bansal, 2016; Li et al., 2017, 2019). Most
of these works have utilized different variants of long short term memory
(LSTM) networks for both entity recognition and relation identification tasks.
Though these networks have shown promising results, they fail to capture best
feature representations.

There have been ground-breaking developments in NLP with the advent of
transformer architectures like BERT (Devlin et al., 2018), GPT (Radford et
al., 2019), RoBERTa (Liu et al., 2019) etc. Transformer networks utilize the
power of multi-head self-attention mechanism to capture context-sensitive
embeddings and interactions between tokens. Since they have been pre-trained
on extremely large documents for language modeling task, they have shown state
of the art results with downstream tasks like sentiment analysis, question
answering, machine translation etc.

Our work exploits the use of ”RoBERTa” architecture for ADE task using a
question-answering(QA) framework. There are three components that lead to the
final identification of ADE: NER module, Classification module to identify
relevant texts and most important component i.e. a Query module fine-tuned for
establishing relation between drug and the ADE. There has been very limited
work in this direction that leverages the transformer architecture in QA setup
to extract entity-relation especially in context of ADE identification using a
drug.

Rest of our work is organized as follows. Section 2 contains related work. In
Section 3, we discuss our approach in detail i.e. system architecture,
dataset, experimental setup, training and evaluation metrics. Section 4
discusses the experimental results which is followed by a conclusion in
Section 5.

##  2\. Related Work

Initial works in NLP focussed on the use of statistical features or kernel-
based methods for relation extraction or classification related tasks (Zelenko
et al., 2003; Chan and Roth, 2011). These research works focus on a pipeline
approach where entities are identified first using a NER module and then a
relationship is identified between those entities. Explicitly computing these
statistical features poses a requirement for a large amount of annotated
documents for training purposes and they are also associated with risk of
errors propagating to different steps.

Li & Ji (Li and Ji, 2014) tackled the joint entity relation extraction problem
as a single transition based model using integer linear programming method to
draw inferences. Deep Learning methods have recently picked up a lot of pace
and contributed a lot in field of NLP. In biomedical domain, the task of NER
to identify drug, disease, dosages etc. has been worked out as a sequence
labelling problem where tokens are tagged according to BILOU/BIO scheme. Miwa
& Bansal (Miwa and Bansal, 2016) proposed an end-to-end joint relation
extraction model where they stacked bidirectional tree-structured LSTMs on
bidirectional sequential LSTMs. Their work inspired many researchers to take
this knowledge to biomedical entity and relation extraction tasks.

Li et. al (Li et al., 2017) proposed a neural joint model for ADE relation
extraction from biomedical texts. They utilized Bi-LSTM architecture for
biomedical entity recognition using concatenation of character representation,
POS embedding and word embedding as input features. This layer shared partial
parameters with another Bi-LSTM layer which was tuned during a joint training
process. These works using Bi-LSTMs helped to make a significant progress in
entity-relation extraction tasks. However, the multi-head self attention
capabilities allowed transformers to capture long range dependencies
efficiently. This along with contextual dense representations from pre-
training using a very large corpus allowed a better feature extraction
capability.

Li et. al (Li et al., 2019) in their work formalize a QA framework using BERT
architecture. Their research emphasizes how question based queries can encode
the important information for entity-relation class identification in the
question and at the same time provide a natural way of jointly modeling entity
and relation. Eberts and Ulges (Eberts and Ulges, 2019) present a span-based
joint entity and relation extraction model with transformer pre-training. They
add a span classification layer that filters entities from non-entities. The
filtered entities are then used for relation identification purposes. We build
upon the learnings of these research endeavours to devise a new QA framework
for ADE identification task using a more powerful transformer architecture
RoBERTa that we describe in the next section.

##  3\. Approach

In this section, we introduce our system architecture (Figure 1) and explain
different modules it invokes.

Figure 1. End-to-End System Architecture (a) Positive (ADE) and Negative (Non-
ADE) sentences are passed through NER module (b) NER module filters sentences
which don’t have drug entity mentions (c) Selected sentences at NER stage are
passed through Classification Module which filters out sentences where Drug-
ADE relation is less probable (d) RoBERTa QA Module identifies ADE
relationship (e.g. depressed level of consciousness) corresponding to the
selected drug entity (e.g. ibuprofen)

###  3.1. Entity Recognition Module

Name entity recognition has been identified as a pivotal task in NLP.
Classification of words from biomedical text into predefined categories like
drug, disease, dosage etc. is a challenging problem. Many researchers have
faced an issue of unavailability of annotated medical corpus due to which
model generalization becomes difficult.

In our system, to identify the drug names in a given phrase, we leveraged
recently developed Med7 (Kormilitzin et al., 2020) NER module which is trained
on a collection of 2 million free-text patients’ record from MIMIC-III corpus
followed by fine-tuning on the NER task. ”It has attained a lenient micro-
average F1 of 0.957 across seven different categories (Dosage, Drug, Duration,
Form, Frequency, Route, Strength)”.

###  3.2. Classification Module

Denoising and extraction of relevant information in textual data is a crucial
step as it improves the learning mechanism and generalizing capability of any
model. After recognizing the drug entities (section 3.1), to identify the
phrases where at least one drug and adverse event pair coexists at stage 2, we
trained a Bi-LSTM (Schuster and Paliwal, 1997) based binary classifier on ADE
sentences (section 4.1) and cross-validated it in K-Fold setting.

This model aims to filter out phrases with no presence of drug and adverse
event pair and helps in improving the performance of our Q&A module and
thereby making it more robust.

Figure 2. Q&A Module for Adverse Drug Event Identification. Input segment A
consists of a drug acting as an alias for question and segment B consists of
context where an ADE is mentioned. Encoded representation of segments passed
as an input through 12-layered transformer network. 1-D convolutional layer is
applied on top of 768 dimensional embedded representations followed by a
softmax activation to identify the ADE.

###  3.3. Q&A Module

Given a passage of text with a user query, a question-answering system
discovers the span of text in the passage that best describes answer to the
question being asked.

BERT (Devlin et al., 2018) based pre-trained language models have achieved
state-of-the-art performances on multiple tasks like semantic role labeling,
question answering, machine translation, etc as it has been trained on large
scale corpora and generalizes the downstream task very well. It can learn the
extremely complex representation in text with a transformer-based (Vaswani et
al., 2017) self-attention mechanism and play a crucial role in improving
varieties of NLP systems.

Liu et. al (Liu et al., 2019) in their paper that introduces RoBERTa, observed
BERT to be ”significantly undertrained” and also found out that with a better
selection of hyperparamerters and training size, its performance could be
considerably improved. RoBERTa, a re-implemented version of BERT in FAIRSEQ
(Ott et al., 2019), has been trained with different learning rate, number of
warm-up steps, batch size and out-performs state-of-the-art results on GLUE,
RACE and SQuAD dataset. We use a pre-trained RoBERTa base model and fine-tune
it on drug-related adverse effects corpus (section 4.1) to identify the
adverse event corresponding to a drug by adding a Q&A head (Figure 2). First
we process the data in desired input format for RoBERTa into 2 segments A and
B. Segment A consists an encoded vector of drug treated as a question followed
by segment B that consists another encoded vector of context/sentence where
adverse event is mentioned. Then we pass this processed data into a 12-layered
transformer network of RoBERTa and use its output that represents the 768
dimensional learnt embeddings of encoded input for further processing. After
this, we apply a one dimensional CNN layer with a (1 x 1) convolution filter
that creates a feature map of these embeddings followed by a softmax
activation layer to predict the probability of start/end tokens of the adverse
event present in a span of the given text.

##  4\. Experimental Studies

###  4.1. Dataset and Evaluation Metrics

We use drug-related adverse effects corpus (Gurulingappa et al., 2012)
containing sentences from 1644 PubMed abstracts. These sentences are divided
into 2 categories (i) ADE (ii) Non-ADE. Former consists of sentences where at
least one pair of drug and its adverse effect is present while latter consists
of sentences with no such pair.

Table 1. Distribution of sentences in dataset Category |  Number of Unique Sentences
---|---
ADE | 6617
Non-ADE | 16688

Examples for ADE and Non-ADE instances are:

  * •

ADE: 14-year-old girl with newly diagnosed sle developed a pruritic bullous
eruption while on prednisone

  * •

Non-ADE: This patient did not have any predisposing factors for the
development of an aortic thrombus before the chemotherapy was initiated.

To gauge the performance of our system, we use common performance metrics such
as Precision, Recall and F1.

(1) $P=\frac{TP}{TP+FP},\\\ R=\frac{TP}{TP+FN},\\\ F1=\frac{2*P*R}{P+R}$

###  4.2. Training

We perform training for 2 different modules involved in our system i.e.
Classification and Q&A.

(i) For classification module, we construct a train and test dataset by
selecting randomly sampled 9931 training and 1272 testing instances. We train
a Bi-LSTM with Adam optimizer (Kingma and Ba, 2014) minimizing binary-
crossentropy loss in stratified K-Fold (k=10) setting to ensure consistency in
our model performance. For final prediction on hold out dataset, we use an
ensemble of these 10 Bi-LSTMs. Below appears the distribution of target
variable in train and test sets:

Table 2. Distribution of target variable in train and test for Classification Module Positive (ADE) |  Negative (Non-ADE ) | Dataset
---|---|---
3976 |  5955 | Train
610 |  662 | Test

(ii) To predict ADE corresponding to a drug, we create the train & and test
data by determining 5955 training and 662 testing instances from ADE sentences
(section 4.1) with random sampling.

We use pre-trained RoBERTa base model and finetune it on these 5955 sentences
with a 1D CNN head (Figure 2) followed by a softmax activation layer to
generate probability corresponding to start and end token of adverse event
present in a context.

We train this entire architecture in K-Fold setting (k=5) with 3 epochs each
on 12GB Nvidia P100 GPU. We use an ensemble of prediction probabilities for
start & end tokens generated by models trained on each of the 5 folds. We set
a learning rate of 3e-5 for Adam optimizer and categorical-crossentropy as our
loss function with a label smoothing of 0.1.

##  5\. Results

Our RoBERTa based QA approach utilizes drug entity passed as a question to
determine answer i.e. adverse drug event in the given context. We determine
the performance at two levels: (i) Performance of individual modules (ii)
Performance of entire architecture (system).

Previous works (Li et al., 2017, 2019; Eberts and Ulges, 2019) in the field
focus on approaches that leverage either a joint training approach or a
cascading pipeline approach to identify entities and then classify those
extracted entities to ascertain existence of any relationship. However, our
approach as described in previous section is not tailored similar to these
settings. Hence, a direct comparison of individual modules is difficult. We
elucidate our approach for calculation of system’s performance using these
three modules together.

###  5.1. Performance of Individual Modules

Table 3 below details the performance of Classification and QA module
corresponding to their selected training and validation set as described in
the previous section. The results for classification module are reported based
on mean results from 10-Fold validation sets. The presence of Non-ADE
sentences (Section 4.1) where drug is either not known or wrongly identified
is expected to create a bias in understanding of effectiveness of RoBERTa QA
module for ADE identification tasks. Hence, we use only ADE sentences where
drug is known beforehand for determining the performance of Roberta QA
framework. In this scenario, the precision would equal to 1 and hence we use
recall to gauge the true efficacy of QA module.

Table 3. Performance Evaluation of Classification & QA Modules Metrics |  Classification Module | QA Module
---|---|---
Precision |  82.74 | -
Recall |  81.44 | 87.37
F1 |  82.06 | -

###  5.2. Performance of End-to-End Architecture

Errors generated by different components in a system create a cascading effect
and this aggregation of errors might render the system to a practically futile
state.

In real-world applications, the overall task comprises of drug identification,
noise removal and then drug-ADE relationship identification. In that process,
obliteration of noisy textual information should also be accounted into the
success criteria for smooth functioning of a NLP pipeline. We describe our
approach to calculate the efficacy of entire system through Figure 3 below.

Figure 3. End-to-End System Architecture’s Overall Performance Calculation.
Sentences in two sets (Spos & Sneg) are passed into the system with an
objective to identify drug-ADE relationship from Spos and removal of all the
Sneg instances. Classification matrix generated at each stage is leveraged for
final calculation of precision, recall and F1 metrics.

At each stage either a sentence is eliminated from system or it moves to next
module for operation. Spos [+D,+C,+A] denotes that among the positive
sentences passed as input, these sentences had drug and were sent to QA module
by classification module where ADE was accurately identified by RoBERTa QA
architecture. Sneg[-D], Sneg[+D,-C] contribute to correctly removed instances
from system. Similar analogy for their Spos counterparts identifies
misclassified samples from the system. Equations for calculation of Precision,
Recall and F1 scores for the entire architecture can be calculated using (1)
and the modified definitions below:

$\displaystyle TP=S_{neg[-D]}+S_{neg[+D,-C]}+S_{pos[+D,+C,+A]}$ (2)
$\displaystyle FN=S_{pos[-D]}+S_{pos[+D,-C]}$ $\displaystyle
FP=S_{pos[+D,+C,-A]}+S_{neg[+D,+C,-A]}$

We detail the final performance metrics for the end-to-end system architecture
in Table 4. We also visit the effectiveness of prominent approaches for ADE
identification task in Table 5 below. After a thorough study of the relevant
literature for ADE identification tasks, we observed that in joint training
approaches, reporting of results is done in a way that the relation
classification task also incorporates effects due to errors generated by
entity recognition modules. Hence, we compare the results of relation
classification task to the effectiveness of our entire end-to-end
architecture.

Table 4. End-to-End System Performance  P |  R | F1
---|---|---
88.37 |  84.44 | 86.36
Table 5. Comparison with different methods Methods |  P |  R |  F1
---|---|---|---
CNN + Global features (Li et al., 2016) |  64.00 |  62.90 |  63.40
BiLSTM + SDP (Li et al., 2017) |  67.50 |  75.80 |  71.40
SpERT (Eberts and Ulges, 2019) |  77.77 |  79.96 |  78.84
Our Model (End-to-End) |  88.37 |  84.44 |  86.36

With a closer and detailed inspection, we observe that even after error
propagation effect, the overall effectiveness of our system is better than
existing approaches for ADE identification.

##  6\. Conclusion

In this paper, we propose a novel approach for ADE identification tasks in
biomedical texts. The experimental results highlight the effectiveness of
using an end-to-end pipeline comprising of NER, classification and RoBERTa
based QA modules where it achieves competitive performances in comparison to
the existing best systems. We hope that our work would open up a new dimension
in ADE identification and entity-relation tasks in general.

Although our end-to-end architecture achieves promising results, we plan to
improve on the NER component in the architecture by building a transformer
based biomedical NER. We also plan to extend this pipeline to social media
platforms like twitter for pharmacovigilance based studies.

## References

  * (1)
  * Bates et al. (1995) David W. Bates, David J. Cullen, Nan Laird, Laura A. Petersen, Stephen D. Small, Deborah Servi, Glenn Laffel, Bobbie J. Sweitzer, Brian F. Shea, Robert Hallisey, Martha Vander Vliet, Roberta Nemeskal, Lucian L. Leape, David Bates, Patricia Hojnowski-Diaz, Stephen Petrycki, Michael Cotugno, Heather Patterson, Mairead Hickey, Sharon Kleefield, Jeffrey Cooper, Ellen Kinneally, Harold J. Demonaco, Margaret Dempsey Clapp, Theresa Gallivan, Jeanette Ives, Kathy Porter, B. Taylor Thompson, J. Richard Hackman, and Amy Edmondson. 1995.  Incidence of Adverse Drug Events and Potential Adverse Drug Events: Implications for Prevention.  _JAMA_ 274, 1 (07 1995), 29–34.  https://doi.org/10.1001/jama.1995.03530010043033
  * Chan and Roth (2011) Yee Seng Chan and Dan Roth. 2011.  Exploiting syntactico-semantic structures for relation extraction. In _Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies_. 551–560.
  * Chiatti et al. (2012) Carlos Chiatti, Silvia Bustacchini, Gianluca Furneri, Lorenzo Mantovani, Marco Cristiani, Clementina Misuraca, and Fabrizia Lattanzio. 2012\.  The economic burden of inappropriate drug prescribing, lack of adherence and compliance, adverse drug events in older people.  _Drug safety_ 35, 1 (2012), 73–87.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ (2018).
  * Eberts and Ulges (2019) Markus Eberts and Adrian Ulges. 2019.  Span-based Joint Entity and Relation Extraction with Transformer Pre-training.  _arXiv preprint arXiv:1909.07755_ (2019).
  * Edlavitch (1988) Stanley A Edlavitch. 1988\.  Adverse drug event reporting: improving the low US reporting rates.  _Archives of internal medicine_ 148, 7 (1988), 1499–1503.
  * Gurulingappa et al. (2012) Harsha Gurulingappa, Abdul Mateen Rajput, Angus Roberts, Juliane Fluck, Martin Hofmann-Apitius, and Luca Toldo. 2012.  Development of a benchmark corpus to support the automatic extraction of drug-related adverse effects from medical case reports.  _Journal of biomedical informatics_ 45, 5 (2012), 885–892.
  * Jagannatha et al. (2019) Abhyuday Jagannatha, Feifan Liu, Weisong Liu, and Hong Yu. 2019.  Overview of the first natural language processing challenge for extracting medication, indication, and adverse drug events from electronic health record notes (MADE 1.0).  _Drug safety_ 42, 1 (2019), 99–111.
  * Johnson and Booman (1996) Jeffery Johnson and Lyle Booman. 1996.  Drug-related morbidity and mortality.  _Journal of Managed Care Pharmacy_ 2, 1 (1996), 39–47.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * Kormilitzin et al. (2020) Andrey Kormilitzin, Nemanja Vaci, Qiang Liu, and Alejo Nevado-Holgado. 2020.  Med7: a transferable clinical natural language processing model for electronic health records.  _arXiv preprint arXiv:2003.01271_ (2020).
  * Li et al. (2017) Fei Li, Meishan Zhang, Guohong Fu, and Donghong Ji. 2017\.  A neural joint model for entity and relation extraction from biomedical text.  _BMC bioinformatics_ 18, 1 (2017), 1–11.
  * Li et al. (2016) Fei Li, Yue Zhang, Meishan Zhang, and Donghong Ji. 2016\.  Joint Models for Extracting Adverse Drug Events from Biomedical Text.. In _IJCAI_ , Vol. 2016. 2838–2844.
  * Li and Ji (2014) Qi Li and Heng Ji. 2014\.  Incremental joint extraction of entity mentions and relations. In _Proceedings of the 52nd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_. 402–412.
  * Li et al. (2019) Xiaoya Li, Fan Yin, Zijun Sun, Xiayu Li, Arianna Yuan, Duo Chai, Mingxin Zhou, and Jiwei Li. 2019\.  Entity-relation extraction as multi-turn question answering.  _arXiv preprint arXiv:1905.05529_ (2019).
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.  _arXiv preprint arXiv:1907.11692_ (2019).
  * Miwa and Bansal (2016) Makoto Miwa and Mohit Bansal. 2016.  End-to-end relation extraction using lstms on sequences and tree structures.  _arXiv preprint arXiv:1601.00770_ (2016).
  * Ott et al. (2019) Myle Ott, Sergey Edunov, Alexei Baevski, Angela Fan, Sam Gross, Nathan Ng, David Grangier, and Michael Auli. 2019\.  fairseq: A fast, extensible toolkit for sequence modeling.  _arXiv preprint arXiv:1904.01038_ (2019).
  * Radford et al. (2019) Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2019\.  Language models are unsupervised multitask learners.  _OpenAI Blog_ 1, 8 (2019), 9.
  * Rocchiccioli et al. (2007) Judith T Rocchiccioli, Julie Sanford, and Bonnie Caplinger. 2007\.  Polymedicine and aging. Enhancing older adult care through advanced practitioners. GNPs and elder care pharmacists can help provide optimal pharmaceutical care.  _Journal of gerontological nursing_ 33, 7 (July 2007), 19—24.  https://doi.org/10.3928/00989134-20070701-04
  * Schuster and Paliwal (1997) Mike Schuster and Kuldip K Paliwal. 1997.  Bidirectional recurrent neural networks.  _IEEE transactions on Signal Processing_ 45, 11 (1997), 2673–2681.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need. In _Advances in neural information processing systems_. 5998–6008.
  * Wu et al. (2012) Chen Wu, Chaim M Bell, and Walter P Wodchis. 2012.  Incidence and economic burden of adverse drug reactions among elderly patients in Ontario Emergency Departments.  _Drug safety_ 35, 9 (2012), 769–781.
  * Zelenko et al. (2003) Dmitry Zelenko, Chinatsu Aone, and Anthony Richardella. 2003\.  Kernel methods for relation extraction.  _Journal of machine learning research_ 3, Feb (2003), 1083–1106.


# On upper bounds for the multi-fold
chromatic numbers of the plane

Jaan Parts

(Kazan, Russia<EMAIL_ADDRESS>
###### Abstract

The multi-fold chromatic number of the plane $\chi_{m}$ is the smallest number
of colors $k$, sufficient to color each point of the Euclidean plane in
exactly $m$ colors, so that for any pair of points at a unit distance from
each other, two corresponding $m$-subsets of $k$-set do not contain any common
color. We consider upper bounds for $m$-fold chromatic numbers of the plane.
Our main result is that for any $m$ the inequality
$\chi_{m}<(1+2/\sqrt{3})^{2}\cdot m+3.501$ holds.

##  1 Background

The classical chromatic number (CN) of a plane is introduced as follows. Each
point of the Euclidean plane is associated with some color $c\in C$. A proper
coloring is one in which any pair of points located at a unit distance from
each other has different colors. The smallest number of colors $k=|C|$ that
provides a proper coloring is called the chromatic number of the plane $\chi$.

Even $2/3$ of a century after the definition of CN, its exact value is still
unknown, and only the range of possible values of $\chi$ is known. During this
time, various extensions and generalizations of the concept of CN have been
invented. Now there are so many different generalizations111 There is nothing
easier than coming up with a new CN: it is enough to modify the coloring
rules, change the object of coloring, or introduce an additional condition for
some of its characteristics. that it is difficult not to get confused in them.
One of the examples is the fractional CN (FCN), as well as the multi-fold (or
$m$-fold) CN (MCN). As with CN, the exact values of MCN and FCN are not known.

In the case of MCN $\chi_{m}$, each point of the plane is associated with a
set of $m$ different colors222 One should distinguish between MCN and the list
CN (LCN). For the latter one each point is associated with its own list of
allowed colors, but receives one single color in the end. (the value of $m$ is
the same for all points), while for any pair of points at a unit distance from
each other, the $m$-sets should not have common colors. Clearly, CN is a
special case of MCN with $m=1$. FCN is defined as
$\chi_{f}=\lim_{m\rightarrow\infty}\frac{\chi_{m}}{m}$. That is, MCN is
intermediate between CN and FCN.

In what follows, we are only interested in the upper bounds for MCN (best
known FCN and CN upper bounds are $\chi_{f}<4.3599$ and $\chi\leq 7$). The
general approach for obtaining such estimates is to create some periodic
tiling, in which the plane is densely covered $m$ times (layers) with tiles of
$k$ different colors, so that the maximum width of tiles does not exceed 1,
and a distance of at least 1 is maintained between tiles of the same color. It
remains only to propose the most suitable shape of the tiles.

This problem was previously considered by Grytczuk et al. [3], where regular
hexagons formed the basis of the tiling. We improve the known upper bounds $k$
of $\chi_{m}$ by considering other tile shapes.

##  2 Tools

### Basic constructions.

We use two basic constructions depending on the value of $m$. In both cases,
the plane is colored with the same layers, shifted relative to each other.
Each layer is covered without gaps by identical tiles, the centers of which
form some lattice.

In the case of small $m$, each of the $m$ layers is covered by hexagonal tiles
of $k$ different colors, and can be divided into identical clusters of $k$
tiles (one tile of each color). The layers are overlapped so that tiles of the
same color from different layers are shifted from each other by the distance
$d$ not less than 1.

In the case of large $m$, each of the $k$ layers contains tiles of only one
color and can be partitioned into identical fragments, consisting of
$k=k_{1}\times k_{2}$ small congruent cells, $m$ of which form one colored
tile (with a width not exceeding 1), and the rest remain uncolored. The shape
of the tiles is chosen so as to increase the proportion of the colored plane.
The layers are overlapped so that each point of the plane is colored in
exactly $m$ different colors.

The main task is to minimize the $k/m$ ratio for each $m$.

### Asymptotic estimates.

The optimal shape of the tiles depends on the number of layers allowed. For
small $m$, the tiles should be angular to ensure a uniform coverage of the
plane without gaps. Hexagonal tiles are best suited for this purpose. The
smallest value of the $k/m$ ratio, which can be achieved theoretically by
tiling the plane with regular hexagons with sides $1/2$, located at a unit
distance from each other at the nodes of a hexagonal lattice, is
$\tau_{1}=(1+2/\sqrt{3})^{2}\approx 4.642734410$. The density of the tiling
depends on the angle $\alpha$ of the rotation of the hexagons relative to the
lattice (or, equivalently, on the shift between the parallel sides of adjacent
tiles of the same color). So the smallest $k/m$ value changes from $\tau_{1}$
to $16/3$ for extreme cases, where tiles are turned to each other by sides or
corners, respectively. More precisely, $k\geq\lceil\tau_{1}\cdot m\rceil$, due
to the integer value of $k$.

At sufficiently large $m$, it becomes possible to increase the plane coloring
density by using more rounded tiles [4]. So,when using circle tiles, $k/m$ can
approach the threshold $\tau_{2}=8\sqrt{3}/\pi\approx 4.410631163$. The best
known bound
$\tau_{3}=\min_{\,\theta}\left[(1+\cos\theta)^{2}/\sqrt{3}\,/\left(\pi/6-\theta+\sin(2\theta)/2\right)\right]\approx
4.359868202$ is given by Croft’s construction [2] (a rounded hexagon333 In
publications, one can also find such names as roundgon, or tortoise. with
equal lengths of the straight and rounded parts of the perimeter).

One can search for $k$ in form $k=\tau_{1}\cdot m+\delta$ or
$k=\lceil\tau_{1}\cdot m\rceil+\Delta$ with some real or integer correction
value $\delta$ or $\Delta=\lfloor\delta\rfloor$ respectively.

### Tile shape.

We will focus basically on the case of small $m$, where hexagonal tiles are
used. Here the problem is to find for each number of layers $m$ a tiling that
minimizes the number of tiles $k$ in the repeating cluster.

To simplify the construction of the tiling, additional restrictions are
imposed on the parameters of the tiles: a) all 6 corners of the tiles are
located on a circle of unit diameter; b) all distances $d$ between the nearest
three tiles of the same color (from different layers) are the same.

Thereby, each tile is symmetrical about the center, has three pairs of
parallel sides and is characterized by only two free parameters, for example,
the width $w$ and the position of the upper corner $u$ (see Fig.1). Looking
ahead, we note that further complication of the tile shape (with one
exception) and the introduction of additional free parameters, in our
experience, did not lead to a decrease in $k$. The parallel sides of the tile
are the same and form a rectangle with unit diagonals. We place one of these
rectangles vertically, so that a horizontal row of tiles is formed.

 |
---|---
Figure 1: To the definition of the tiling parameters by the example of
$m=3\times 2=6$ and $k=31$. The base layer tiles are shown, as well as the
base color tiles (highlighted in gray). Other layers are hidden so as not to
clutter up the picture, but can be unambiguously restored by the position of
the base color tiles. A possible clustering option is shown in bold lines. The
enlarged fragment of the tiling on the right shows the three nearest tiles of
the same color on a square grid with a step 0.1.

### Tiling parameters.

Some color (tile, cluster, layer) is declared base.

For the convenience of indexing and determining the position of the tiles,
oblique coordinates $p_{n}=(i_{n},j_{n})$ are introduced, where $i_{n}$ and
$j_{n}$ are integer indices of the tile $n$ along two (out of three possible)
directions of the hexagonal lattice coinciding with the centers of the tiles
in the base layer. The tiles of the remaining layers generally have fractional
indices.

Adjacent tiles of the base color (from different layers) form another
hexagonal lattice and are shown in gray in the figures. More specifically, $m$
adjacent tiles (one from each layer) form a parallelogram corresponding to one
of the possible variants of integer decomposition $m=r\times s$.

To carry out the calculations, it is enough to consider the coordinates of six
tiles of the base color, numbered from 0 to 5 (see Fig.1): $p_{0}=(0,0)$ for
the base tile; $p_{1}$ and $p_{2}$ for the base tiles of two adjacent clusters
of the same layer; $p_{1}$ and $p_{3}$ for the tiles corresponding to the
decomposition $m=r\times s$; $p_{0}$, $p_{4}$ and $p_{5}$ for adjacent tiles
that form a triangle of the minimum perimeter; moreover, $p_{5}=p_{1}/r$,
$p_{4}=p_{3}/s$, $p_{3}=p_{2}-t\cdot p_{5}$.

The concrete construction of the tiling is described by the parameters
characterizing the density of the tiling $\\{m,k\\}$, the tile shape
$\\{w,u\\}$, the distance between the tiles $\\{d\\}$, the shape of the
cluster $\\{i_{1},j_{1},i_{2},j_{2}\\}$, an arrangement of tiles of the same
color $\\{r,s,t\\}$. Note that these parameters are not fully independent.

### Computational procedure.

We seek to minimize $k$ and, subject to that, to maximize $d$ for each $m$.
The search for the optimal tiling is performed by exhaustive search. For given
$m$ the values of $k$ are sequentially increased, starting from
$k=\lceil(1+2/\sqrt{3})^{2}\cdot m\rceil$, until a tiling is found with $d\geq
1$. For each $m$, all possible variants of $\\{r,s,t\\}$ are checked. For each
$k$, a set of possible $\\{p_{1},p_{2}\\}$ is determined taking into account
the relation $k=i_{1}\cdot j_{2}-i_{2}\cdot j_{1}$. For each configuration
described by integer parameters $\\{m,k,i_{1},j_{1},i_{2},j_{2},r,s,t\\}$, a
procedure is started with variable parameters $w$ and $u$, which maximizes the
value of $d$, so that $0\leq u\leq w\leq 1\leq d$.

To reduce the search area, a number of other natural restrictions are
introduced on the range of parameter values: $j_{1}\leq 0$, $i_{1}\geq
2\,|j_{1}|$ (take into account the symmetry of the tiling); $i_{1}>2\,r$
(provides a horizontal gap between the tiles); $k>2\,s\,i_{1}$ (sets a lower
limit on the cluster size); $0\leq i_{2}<i_{1}$, $j_{2}\leq 2\,s+1$,
$i_{2}+j_{2}>2\,s$, $t<r$ (limit the position $p_{2}$).

##  3 Tricks

### Extrapolation.

If one superimpose two arbitrary tilings with parameters $\\{m_{1},k_{1}\\}$
and $\\{m_{2},k_{2}\\}$ one on top of the other, then one can always get a
combined tiling with parameters $\\{m,k\\}$ using $k=k_{1}+k_{2}$ different
colors and $m=m_{1}+m_{2}$ layers.

Combining the tilings allows one to obtain improved MCN estimates for some $m$
from the available MCN values for smaller $m$, as well as to extrapolate to
larger $m$ values. For this, all variants of the decomposition of $m$ into
terms $m=m_{1}+m_{2}$ are sorted out, and for them the minimum sum
$k=k_{1}+k_{2}$ is determined. Practice shows that such extrapolation gives
good results (see Section 4).

### Wavy sides.

If the parallel sides of adjacent tiles of the same color are shifted by more
than half the length of a side444 Notice that if all tiles have the same shape
then each side of the tile must be symmetrical about its center. Thus the
position of the center relative to the ends of the side is fixed., one can use
the trick of drawing the relevant part of this side as an arc with a radius of
1 centered at the nearest corner of the adjacent tile to slightly reduce the
distance between the tiles. The largest effect (reduction of distance by about
$0.8\%$) is observed when the sides of adjacent tiles are shifted by $3/4$ of
their length.

Potentially, the effect of wavy sides can work independently for any pair of
parallel sides. Since, with increasing $m$, there is a general tendency
towards a decrease in the angle $\alpha$ (tiles tend to take the shape of a
regular hexagon and turn sides to each other), one can expect the appearance
of wavy sides only for small $m$ (practically for $m<40$).

##  4 Results

### Small $m$.

Let’s start with the first construction. Here for a particular pair
$\\{m,k\\}$, it is sufficient to get the only tiling with $d\geq 1$. Usually,
one can find several different tilings with this property555 Also note that
the chosen notation system has redundancy, and the same tiling can be
specified in different ways.. We show those that maximize $d$. Along with $k$,
we use the deviation $\delta$ from the theoretical minimum $\tau_{1}$ for
hexagonal tiles.

We did a full search for all $m\leq 100$, as well as for $m\leq 2000$ with
some simplifications and gaps. Further search is devoid of practical sense,
bearing in mind the increase in its labor intensity, the possibility of simple
extrapolation, as well as the transition to the second construction.

Most of the obtained values $\delta=k-\tau_{1}$ are within the range
$\delta\in(0.35,3)$. The exceptions are $m\in\\{7,48,132,135\\}$. For all
$m>135$ we got $\delta<2.74$. The largest value $\delta<3.501$ is observed at
$m=7$.

The lower threshold $\delta>0.35$ is explained by the fact that, if the shape
of the tiles approaches a regular hexagon, a row of tiles of the same color
($r$-axis on Fig.1) should have a slight slope with respect to a row of tiles
of the same layer ($i$-axis). This means that the value $\tau_{1}$ is
unattainable in this construction for finite $m$.

Oddly enough, we managed to find the only (but practically important) case
when the use of wavy sides leads to a decrease in $k$: for $m=6$, it is
possible to reduce $k$ from 31 to 29 (due to reduction of $d$ by only
$0.15\%$).

Figures 2 and 3 show some examples of tilings that give the best found upper
bounds of MCN for $m\leq 10$. If one look closely, one can see the wavy sides
in the case $k/m=29/6$. The discovery of this tiling increases the efficiency
of the extrapolation trick. In the range $100<m\leq 2000$, we got 1456
explicit tilings with $\delta<3$. The remaining 444 estimates of $\delta$ are
obtained by extrapolation. Besides, extrapolation allowed to decrease $\delta$
by 1 and 2 for 87 and 12 cases of explicit tiling, respectively.

7/1 |   |
---|---|---
11/2 |   |
16/3 |   |
21/4 |   |
26/5 |   |
Figure 2: Tilings for $m=1\ldots 5$ and best found values of $k/m$. 29/6 |   |
---|---|---
36/7 |   |
39/8 |   |
43/9 |   |
49/10 |   |
Figure 3: Tilings for $m=6\ldots 10$ and best found values of $k/m$. Table 1: Characteristics of the best found tilings. $m$ |  $k$ |  $k^{*}$ |  $k/m$ |  $\delta$ |  $i_{1}$ |  $j_{1}$ |  $i_{2}$ | $j_{2}$ |  $r$ |  $s$ | $t$ |  $w$ |  $u$ |  $d$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 |  7 |  7 |  7.0000 |  2.3573 |  3 |  -1 |  1 | 2 |  1 |  1 | 0 |  0.866025 |  0.433013 |  1.3229
2 |  11 |  12 |  5.5000 |  1.7145 |  5 |  -1 |  1 | 2 |  2 |  1 | 0 |  0.780412 |  0.231291 |  1.0550
3 |  16 |  16 |  5.3333 |  2.0718 |  8 |  0 |  3 | 2 |  3 |  1 | 1 |  0.607271 |  0.303636 |  1.0121
4 |  21 |  24 |  5.2500 |  2.4291 |  5 |  -1 |  1 | 4 |  2 |  2 | 0 |  0.866025 |  0.433013 |  1.0825
5 |  26 |  32 |  5.2000 |  2.7863 |  13 |  0 |  3 | 2 |  5 |  1 | 1 |  0.626508 |  0.272476 |  1.0024
6 |  29 |  32 |  4.8333 |  1.1436 |  7 |  -1 |  1 | 4 |  3 |  2 | 0 |  0.831425 |  0.322944 |  1.0009
6 |  29 |  32 |  4.8333 |  1.1436 |  7 |  -1 |  1 | 4 |  3 |  2 | 0 |  0.831379 |  0.327403 |  0.9994
6 |  31 |  32 |  5.1667 |  3.1436 |  7 |  -1 |  3 | 4 |  3 |  2 | 1 |  0.894452 |  0.362876 |  1.0716
7 |  36 |  37 |  5.1429 |  3.5009 |  16 |  -3 |  12 | 0 |  7 |  1 | 5 |  0.871616 |  0.243410 |  1.0163
8 |  39 |  45 |  4.8750 |  1.8581 |  19 |  -2 |  10 | 1 |  8 |  1 | 4 |  0.792258 |  0.321035 |  1.0091
9 |  43 |  49 |  4.7778 |  1.2154 |  7 |  -1 |  1 | 6 |  3 |  3 | 0 |  0.866025 |  0.433013 |  1.0104
10 |  49 |  55 |  4.9000 |  2.5727 |  12 |  -1 |  1 | 4 |  5 |  2 | 0 |  0.765621 |  0.310729 |  1.0097
11 |  53 |  72 |  4.8182 |  1.9299 |  25 |  -4 |  7 | 1 |  11 |  1 | 3 |  0.928440 |  0.477605 |  1.0080
12 |  57 |  63 |  4.7500 |  1.2872 |  9 |  -1 |  3 | 6 |  4 |  3 | 1 |  0.885280 |  0.397978 |  1.0071
13 |  63 |  87 |  4.8462 |  2.6445 |  30 |  -3 |  21 | 0 |  13 |  1 | 9 |  0.866025 |  0.433013 |  1.0326
14 |  67 |  80 |  4.7857 |  2.0017 |  31 |  -4 |  9 | 1 |  14 |  1 | 4 |  0.931056 |  0.445051 |  1.0034
15 |  71 |  77 |  4.7333 |  1.3590 |  11 |  -1 |  5 | 6 |  5 |  3 | 2 |  0.898147 |  0.379766 |  1.0018
16 |  77 |  81 |  4.8125 |  2.7162 |  9 |  -1 |  5 | 8 |  4 |  4 | 2 |  0.911429 |  0.475233 |  1.0205
17 |  81 |  111 |  4.7647 |  2.0735 |  39 |  -4 |  30 | -1 |  17 |  1 | 13 |  0.872477 |  0.475928 |  1.0171
18 |  86 |  91 |  4.7778 |  2.4308 |  14 |  -2 |  1 | 6 |  6 |  3 | 0 |  0.866025 |  0.433013 |  1.0104
19 |  91 |  123 |  4.7895 |  2.7880 |  43 |  -3 |  16 | 1 |  19 |  1 | 7 |  0.866025 |  0.433013 |  1.0256
20 |  95 |  99 |  4.7500 |  2.1453 |  22 |  -1 |  7 | 4 |  10 |  2 | 3 |  0.868899 |  0.341116 |  1.0086
21 |  99 |  112 |  4.7143 |  1.5026 |  48 |  -5 |  39 | -2 |  21 |  1 | 17 |  0.875384 |  0.503138 |  1.0057
22 |  103 |  120 |  4.6818 |  0.8598 |  49 |  -3 |  18 | 1 |  22 |  1 | 8 |  0.861306 |  0.398575 |  1.0027
23 |  109 |  150 |  4.7391 |  2.2171 |  51 |  -2 |  29 | 1 |  23 |  1 | 13 |  0.858033 |  0.367470 |  1.0126
24 |  113 |  117 |  4.7083 |  1.5744 |  53 |  -3 |  20 | 1 |  24 |  1 | 9 |  0.878789 |  0.419435 |  1.0094
25 |  119 |  121 |  4.7600 |  2.9316 |  11 |  -1 |  9 | 10 |  5 |  5 | 4 |  0.922862 |  0.495630 |  1.0083
26 |  123 |  144 |  4.7308 |  2.2889 |  29 |  -2 |  18 | 3 |  13 |  2 | 8 |  0.883319 |  0.478525 |  1.0135
27 |  127 |  140 |  4.7037 |  1.6462 |  20 |  -1 |  7 | 6 |  9 |  3 | 3 |  0.866025 |  0.433013 |  1.0104
28 |  131 |  144 |  4.6786 |  1.0034 |  61 |  -2 |  35 | 1 |  28 |  1 | 16 |  0.871101 |  0.367408 |  1.0004
29 |  137 |  189 |  4.7241 |  2.3607 |  65 |  -4 |  18 | 1 |  29 |  1 | 8 |  0.869685 |  0.484585 |  1.0128
30 |  141 |  143 |  4.7000 |  1.7180 |  33 |  -1 |  9 | 4 |  15 |  2 | 4 |  0.864112 |  0.407710 |  1.0098
31 |  145 |  201 |  4.6774 |  1.0752 |  69 |  -5 |  29 | 0 |  31 |  1 | 13 |  0.882489 |  0.495857 |  1.0018
32 |  151 |  162 |  4.7188 |  2.4325 |  70 |  -1 |  11 | 2 |  32 |  1 | 5 |  0.864953 |  0.409509 |  1.0143
33 |  155 |  168 |  4.6970 |  1.7898 |  72 |  -1 |  11 | 2 |  33 |  1 | 5 |  0.863131 |  0.387042 |  1.0083
34 |  159 |  184 |  4.6765 |  1.1470 |  75 |  -2 |  42 | 1 |  34 |  1 | 19 |  0.854345 |  0.418193 |  1.0056
35 |  165 |  176 |  4.7143 |  2.5043 |  76 |  -1 |  13 | 2 |  35 |  1 | 6 |  0.876328 |  0.469624 |  1.0131
36 |  169 |  169 |  4.6944 |  1.8616 |  13 |  0 |  0 | 13 |  6 |  6 | 0 |  0.866025 |  0.433013 |  1.0104
37 |  173 |  240 |  4.6757 |  1.2188 |  80 |  -1 |  13 | 2 |  37 |  1 | 6 |  0.875603 |  0.426892 |  1.0061
38 |  177 |  205 |  4.6579 |  0.5761 |  82 |  -1 |  13 | 2 |  38 |  1 | 6 |  0.874272 |  0.406629 |  1.0016
39 |  183 |  203 |  4.6923 |  1.9334 |  83 |  -1 |  17 | 2 |  39 |  1 | 8 |  0.901827 |  0.448498 |  1.0059
40 |  187 |  198 |  4.6750 |  1.2906 |  85 |  -1 |  17 | 2 |  40 |  1 | 8 |  0.900383 |  0.428203 |  1.0022

The best results for $m\leq 40$ are listed in Table 1. Here, for each $m$, the
best found values of $k$ and $d$ are shown, as well as the corresponding
tiling parameters $\\{i_{1},j_{1},i_{2},j_{2},r,s,t,w,u\\}$. For comparison,
the estimates of $k^{*}$ are given according to the data from [3]. They are
not only weaker, but also do not provide convergence to some threshold: for
arbitrarily large $m$ there are estimates of $k^{*}/m$ here both approaching
$\tau_{1}$ (if $m$ has large divisors) and exceeding $3+2\sqrt{3}\approx 6.46$
(if $m$ is prime). For $m=6$, three solutions are shown: one with wavy sides
and two with straight sides with $k=29$ and $k=31$ (option $k=30$ does not
give valid solutions even with wavy sides). The second solution with $k=29$ is
not valid, since it does not satisfy the condition $d\geq 1$, and is presented
for comparison only.

---
Figure 4: The best found values $\delta=k-\tau_{1}\cdot m$ for different $m$
(for two different $m$ scales). Points show the actual values of $\delta$,
circles – their probable estimates $\hat{\delta}$ (see Section 5). The color
indicates that the result belongs to a certain class. The sequences $m\bmod
3=\\{0,1,2\\}$ are shown in black, red and gray, respectively. Prime $m$ are
highlighted in blue. Yellow and green show the values obtained and improved,
respectively, by tricks (wavy sides and extrapolation). Colors are listed in
ascending order of priority. The orange lines correspond to
$\delta=\\{1/3,8/3\\}$.

### Medium $m$.

The results for $m\leq 150$ and $m\leq 1040$ are shown in Fig.4 in the form of
dependence $\delta(m)$. Here the periodicity in $m$ is clearly observed. One
can notice periods of about 3 and 41.75. The both are explained by the fact
that the threshold $\tau_{1}$ is close to $14/3$, so that a) an increase in
$m$ by 3 does little to change the fractional part $f$ of the $k/m$ ratio, b)
$f$ has period $1/(14/3-\tau_{1})\approx 41.78461$.

With further extrapolation of more or less known $k$ for all $m\leq 2000$, it
is possible to obtain estimates $\Delta=\lfloor\delta\rfloor\leq 2$ for
$m<5019$ and $\Delta\leq 3$ for $m<11829$. If we restrict the known range to
$m\leq 500$, then we get $\Delta\leq 2$ for $m<1425$ and $\Delta\leq 3$ for
$m<2704$. If we initially take $m\leq 40$, then after extrapolation
$\Delta\leq 4$ for $m<101$ and $\Delta\leq 15$ for $m<893$. In Fig.4, the
values obtained or improved by extrapolation are highlighted in color.

Figure 5: The tiling with $m=860$ and $k/m<\tau_{1}$.

### Large $m$.

In the case of large $m$ (about 1000 or more), the hexagonal tiles become
ineffective.

For any given threshold value $\tau>\tau_{3}$, one can find $m$ for which the
tiling of the second type will certainly provide $k/m<\tau$. However, such
estimate is rather rough, and for $\tau=\tau_{1}$ gives a value $m$ of the
order $10^{4}$. Alternatively, one can take the densest coverage of the plane
with circles or Croft’s smoothed hexagons, and consider different ways of
filling these tiles with small rectangles (or hexagons), choosing the option
with minimal loss of area around the tile perimeter.

So, for $m=860$, we found a tiling with $k=70\times 57$ (the corresponding
sides of the rectangular cells are $0.0277851\times 0.0306606$) and
$k/m=3990/860\approx 4.63953$ (see Fig.5), for $m=992$ we obtained a tiling
with $k=74\times 62$ and $k/m=4588/992\approx 4.625$. In other words, for
$m=860$ we get $\delta<0$. Combining now the constructions of the first and
second types666 Note that a construction of the second type exists for every
$m$, since, if necessary, one can leave several rectangular cells uncolored.
and using the trick with extrapolation, it is easy to make sure that for all
subsequent $m$ the relation $k/m<\tau_{1}$ also holds, which means that the
previously obtained estimate $\delta\leq 3.501$ is valid for all $m$.

##  5 Puzzles

If we take a closer look at Fig.4, we will notice the following curious thing.
Most of the values of $\delta$ fit into the range $\delta\in[0.35,2.75]$,
moreover, if we divide all $m$ into three equivalence classes ($m\bmod 3$),
then most of $\delta$ fall into the following ranges:

$\delta\in\begin{cases}[1.05,2.05],&m\equiv 0\pmod{3}\\\
[0.75,1.4]\cup[2.4,2.75],&m\equiv 1\pmod{3}\\\
[0.35,0.7]\cup[1.7,2.35],&m\equiv 2\pmod{3}\end{cases}$

If we round off these values to $1/3$, we get the following rough empirical
formula for the probable value of $\delta$:

$\hat{\delta}=f+\begin{cases}1,&m\equiv 0\pmod{3}\\\ \lfloor 3f+m\rfloor\bmod
3,&m\not\equiv 0\pmod{3}\end{cases}$

Here $f=\delta-\lfloor\delta\rfloor=\lceil\tau_{1}\cdot m\rceil-\tau_{1}\cdot
m$ is a fraction part of $\delta$.

For $m\leq 2000$, this rule is violated in about $13.5\%$ of cases. If we
discard the initial and the rough search zones, then in a narrower range
$m\in[150,900]$ the rule is violated in $10.5\%$ of cases. Interestingly,
there are not only positive, but also negative deviations from the estimate
$\hat{\delta}$ (about $3.3\%$ of cases). The question arises: what explains
this rule? First of all, what causes the ”unit holes” in $\delta$ for
$m\not\equiv 0\pmod{3}$? Maybe for such combinations of $k$ and $m$, the tiles
are lined up so badly that for any options of the tiling parameters, they do
not allow to obtain a valid solution. However, such an explanation explains
little, and only states a fact. Something better is required here.

Among the riddles is the only successful case of a trick with wavy sides. One
can try to find other ways to optimize the tiling. For example, is less
regular tiling beneficial with tiles of different shapes in the same layer? Is
it possible to find an intermediate construction for medium $m$, combining an
$m$-layer structure (as for hexagonal tiles) and the tiles of a more complex
shape (like smoothed hexagons), and to cross the threshold $\tau_{1}$ much
sooner, say, closer to $m\approx 100$?

One can notice that the first construction for small $m$ gives convergence to
the threshold $\tau_{1}$ with the rate $1/m$. The second construction for
large $m$ gives a convergence to the smaller threshold $\tau_{3}$, but with
slower rate $1/\sqrt{m}$. Is it possible to improve the convergence in the
second case as well?

The article is finished. But I don’t like to leave the empty spaces.
Therefore, a short story.

Initially, when I had only a few hand-made tilings in my hands, I conceived an
ideal mathematical article with a minimum of text, or maybe with no text at
all to break the record set by Conway and Soifer in their monumental work [1].
In particular, it would immediately relieve me of the need to translate it.
But then, as the material began to accumulate, because my computer also wanted
to take part in the creative process, I had to revise my initial views, and
add some clarifications. Hope you found them interesting. This was facilitated
by A. de Grey and A. Soifer.

## References

  * [1] J.H. Conway and A. Soifer, Can $n^{2}+1$ unit equilateral triangles cover an equilateral triangle of side $>n$, say $n+\varepsilon$? American Mathematical Monthly, no. 112, 2005, p. 78.
  * [2] H.T. Croft, Incidence incidents, Eureka (Cambridge), no. 30, 1967, pp. 22–26.
  * [3] J. Grytczuk, K. Junosza-Szaniawski, J. Sokol, and K. Wesek, Fractional and j-fold coloring of the plane, Discrete & Computational Geometry, vol. 55, no. 30, 2016, pp. 594–609.
  * [4] D. Pritikin, All unit-distance graphs of order 6197 are 6-colorable, Journal of Combinatorial Theory, Series B 73, 1998, pp. 159–163.