# Fine-tuning Transformer-based Encoder for Turkish Language Understanding
Tasks

Savaş Yıldırım<EMAIL_ADDRESS>0000-0002-7764-2891 Istanbul Bilgi
University, Faculty of Engineering and Natural SciencesEski Silahtarağa
Elektrik Santralı Kazım Karabekir Cad. No: 2/13IstanbulTurkey34060

###### Abstract.

Deep learning-based and lately Transformer-based language models have been
dominating the studies of natural language processing in the last years.
Thanks to their accurate and fast fine-tuning characteristics, they have
outperformed traditional machine learning-based approaches and achieved state-
of-the-art results for many challenging natural language understanding (NLU)
problems. Recent studies showed that the Transformer-based models such as
BERT, which is Bidirectional Encoder Representations from Transformers, have
reached impressive achievements on many tasks. Moreover, thanks to their
transfer learning capacity, these architectures allow us to transfer pre-built
models and fine-tune them to specific NLU tasks such as question answering. In
this study, we provide a Transformer-based model and a baseline benchmark for
the Turkish Language. We successfully fine-tuned a Turkish BERT model, namely
BERTurk that is trained with base settings, to many downstream tasks and
evaluated with a the Turkish Benchmark dataset. We showed that our studies
significantly outperformed other existing baseline approaches for Named-Entity
Recognition, Sentiment Analysis, Question Answering and Text Classification in
Turkish Language. We publicly released these four fine-tuned models and
resources in reproducibility and with the view of supporting other Turkish
researchers and applications.

deep learning, transfer learning, transformers, language model

††ccs: Computing methodologies Machine learning algorithms††ccs: Computing
methodologies Natural language processing††ccs: Computing methodologies
Information extraction

##  1\. Introduction

The advent of deep learning and transfer learning in the field of natural
language processing (NLP) has achieved significant improvements for many
tasks. Sequence transduction models such as RNN and LSTM are based on complex
recurrent neural networks. They successfully connected the encoder and decoder
for sequence to sequence problem (Zhang et al., 2020). However one of the
computational bottlenecks of RNN-based models is processing very long
sequences and they have difficulty in keeping long-distance dependencies.

The Transformers (Vaswani et al., 2017) overcome this limitation and simplify
encoder-decoder network architecture by employing attention mechanisms,
dropping recurrence or convolutions entirely. The studies showed that the
Transformer could be trained significantly faster than other standard RNN
architectures since that RNN models are based on more complex recurrent
layers. The initial experiments on some machine translation tasks indicated
that the Transformer-based models are superior in quality. They are more
parallelizable and further requiring significantly less time to train.
Finally, the Transformer-based models became a new paradigm in NLP and led to
many successful derivatives that provides state-of-the-art general-purpose
architectures tested on General Language Understanding Evaluation (GLUE)
Benchmark tasks (Wang et al., 2018); BERT (Devlin et al., 2019), GPT-2
(Radford et al., 2019), RoBERTa (Liu et al., 2019), XLM (Lample and Conneau,
2019), XLNet (Yang et al., 2019), T5 (Raffel et al., 2019) and so forth.

Main advantage of the transformer models is to provide fast and cheap fine-
tuning in terms of space and time complexity by transferring the base language
models trained in the pre-training phase. While the pre-training phase demands
huge amounts of data, the fine-tuning phase requires relatively small amounts
of data depending on downstream tasks. The pre-training phase trains the
language model so that it can build internal representations such as
contextual word embeddings or sentence encoding (Peters et al., 2018), which
allows there after a model to be reused for different downstream tasks even
with little amount of labeled data. The contextualized word embeddings (Peters
et al., 2018), addressed the well known word sense problem and represent the
word within its context contrary to other traditional word embeddings like
word2vec (Mikolov et al., 2013), GloVe (Pennington et al., 2014),
fastText(Bojanowski et al., 2016). Another advantage of the attention
mechanism allows the model to better resolve long-term dependencies and
coreference relation between entities (Joshi et al., 2019). This property is
important especially for machine translation and question answering problems.

In this study we provide a benchmark for Natural Language Understanding (NLU)
downstreams tasks for the Turkish Language. We fine-tuned a pre-built BERT
language model, BERTurk (Schweter, 2020), to four downstream tasks; Named-
Entity Recognition, Sentiment Analysis, Question Answering and Text
Classification. Compared to the baseline approaches, we have got successful
results in the NLU tasks. In order to reproduce the result and support the
community, we publicly released these four fine-tuned models. This study is
considered as the first successful attempt to fine-tune the BERT model for
these four tasks in the Turkish Language.

##  2\. Related Works for Language Models

The traditional language models have been based on n-grams and they were
simply a probability distribution over word sequences with an objective
function of predicting the last n-th word from previous (n-1) words. One
drawback of n-gram language model is that as the model is being trained with a
very large corpus, the number of unique words increases and therefore sparse
data problem arises (Zhang et al., 2020).

Modern neural networks sorted the dimensionality problem out by using
continuous and dense representation for the words, called word embeddings. The
early neural language representation employed encoders to produce static word
embeddings. The most popular models, such as word2vec (Mikolov et al., 2013),
glove (Pennington et al., 2014) and fastText (Bojanowski et al., 2016)
transformed unsupervised textual data into a supervised one by either
predicting the target word using context or predicting contextual words by the
target word based on a sliding window of n-gram context. The final output
turns out to be context-independent word embeddings which ignores word sense.

The main problem of these early models is that they never use the sentence
level representation. Second problem is that the senses of the words as one of
the oldest problems in computational linguistics are ignored and a single
fixed representation is assigned to each word or token. Further studies used a
neural network component which encodes the complete sentence-level information
and contextual word embeddings which finally alleviate word sense problems. In
early days, the essential architectures such as RNN, CNN, LSTM have been
widely used as encoder and decoder in a sequence to sequence problems.

ELMo (Embeddings from Language Models) representations (Peters et al., 2018),
which is based on LSTM, successfully built deep contextualized word
representations and solved the sense and sentence encoding problem. The ELMo
representations differ from traditional word embeddings in that each word
representation is a function of the entire input sentence.

Later on, the Transformers (Vaswani et al., 2017) simplified sequence-to-
sequence models by removing recurrency and only keeping attention. The
Transformer-based architecture such as GPT (Radford et al., 2019) and BERT
(Devlin et al., 2019) successfully captured the contextualized sentence-level
language representations and word embeddings in a modern way. These sentence
encoders differ from sliding window approaches in that they read a complete
sentence instead of any fixed sized word sequences defined by a sliding
window. While Elmo used bidirectional LM to read the entire sentence
independently in both directions and concatenates trained forward and backward
LSTMs, BERT adopted Transformer to read and process the entire input at once
in a simultaneous bidirectional way.

BERT algorithm is able to read the input in a bidirectional way rather than
unidirectional. The BERT model trains two different objectives simultaneously;
Masked Language Model (MLM) and  Next Sentence Prediction (NSP). MLM is based
on optimizing the prediction of the masked words in sentences that are
randomly chosen. MLM randomly replaces some masked tokens with a special
symbol, such as [MASK]. Then, the model training objective is to predict the
masked symbol based on the non-masked context. The major difference between
Transformer-based BERT and LSMT-based ELMo is that the predictions in BERT are
based on the entire context rather than only one direction. The BERT is
however criticized that the masks could be never seen at fine-tuning, which
faces a mismatch between pre-training and fine-tuning phases. NSP is based on
a model that detects whether two sentences follow each other or not. Some
other transformer models such as RoBERTa (Liu et al., 2019) or Albert (Lan et
al., 2020) employ different approaches for the objective functions to improve
the results.

Although both RNN and Transformer architectures gained spectacular
achievements, the main limitation is said to be understanding long-term
dependencies. Transformer-XL (Dai et al., 2019) put together the advantages of
these two architectures. It employs the self-attention modules on each
contextual segment of input and a RNN mechanism to learn dependencies between
successive segments. Based on the inspirations of impressive RNN and
Transformer-based models like ELMo (Peters et al., 2018) and BERT (Devlin et
al., 2019), recently a wide range of LM derivatives have been proposed. The
variants such as  Electra (Clark et al., 2020), Albert (Lan et al., 2020), T5
(Raffel et al., 2019), RoBERTa (Liu et al., 2019) are gaining successful
achievements on GLUE tasks (Wang et al., 2018) which is a set of tools for
evaluating the performance of models across a diverse collection of common NLU
tasks.

##  3\. The Transformer-based Model for the Turkish Language

The advantage of the Transformers is that the same architecture could be used
in both pre-training and fine-tuning except output layers. The pre-trained
model parameters are used to initialize models for different downstream tasks.
Pre-training and fine-tuning made it possible to apply transfer learning in
NLP and have become a new paradigm. In this section, we describe how to fine-
tune BERT-based models for a variety of problems in Turkish; Question
Answering, Sentiment Analysis, Named-Entity Recognition and Text
classification. We also discuss how to evaluate them with a set of benchmark
datasets.

###  3.1. Fine-tuning BERT model for Turkish

We fine-tuned four NLU downstream tasks based on BERTurk (Schweter, 2020) and
uploaded to the repository. BERTurk (Schweter, 2020) has been successfully
pre-trained as a Transformer model with a huge Turkish corpus and shared with
the NLP community. There are two alternatives for BERT training; base and
large. While base architecture has (L=12, H=768, A=12, Total Parameters=110M),
large one has (L=24, H=1024, A=16, Total Parameters=340M) where L is number of
layers, and H is hidden size, A is number of self attention heads. The model
can be used cased or uncased as a way of formatting input text. For the
Turkish language, we have currently only base pre-trained LM of BERT with
cased or uncased.

The last version of Turkish BERT model, which we used in this study, has been
trained on a huge training corpus with file size of 35GB and 45B tokens. The
compiled corpus consists of online resources such as OSCAR Corpus
111https://traces1.inria.fr/oscar/, OPUS corpus 222http://opus.nlpl.eu/,
Wikipedia dump and Turkish Tree-bank 333https://ii.metu.edu.tr/metu-corpora-
research-group. The base model has been trained within a cloud environment
provided by Google TensorFlow Research Cloud (TFRC). Uncased and cased models
with 128K and 32K Vocabulary have been successfully trained. In order to load
or fine-tune the model, we used python and its transformer library
444https://huggingface.co/transformers/, which is built for the NLP community
to share, load and extend their models.

The Transformers can learn syntactic segmentation at the lower levels such as
subwords and then semantic knowledge at the higher levels such as coreference
resolution as they work on specific language. BERT model applies WordPiece
tokenization procedure to each token to understand syntactic knowledge. The
implementation is simply based on the one from tensor2tensor, WordPiece
embeddings (Wu et al., 2016). The wordpiece model is generated using a data-
driven approach to maximize the language-model likelihood of corpus. When a
training corpus is provided, the optimization problem is to select limited
wordpieces. Turkish is an agglutinative language and the root word receives
many suffix so that incredibly large word could be generated as follows
(wikipedia.org)

    "muvaffakiyetsizletiriciletiriveremeyebileceklerimizdenmisiniz"

This word is one of the longest words and can be extended more. It is derived
from the noun word ”muvaffakiyet” (success). It means ”As though you are from
those whom we may not be able to easily make into a maker of unsuccessful
ones”. The response of wordpiece tokenization procedure to the word is truly
successful as follows

    [’muvaffak’,  ’##iyet’,  ’##siz’,  ’##le’, ’##tir’,  ’##ici’,
     ’##le’,  ’##tir’,  ’##iver’,  ’##emeye’,  ’##bilecekleri’, ’##mi’, ’##z’,
     ’##den’, ’##mi’, ’##siniz’, ’##cesine’,  ’##siniz’]

It is now transformed into 18 chunks and represented by (18+2, 32000) tensors
in the transformer systems where +2 is needed for [CLS] and [SEP] tokens. Such
capacity of the syntactic breakdown helps many models to easily fine-tune
downstream tasks. We fine-tuned four models that are now easily accessible via
Hugging-Face platform 555https://huggingface.co/savasy.

###  3.2. Sentiment Analysis

The fine-tuned model has been already shared
666https://huggingface.co/savasy/bert-base-turkish-sentiment-cased. The
dataset for sentiment analysis is obtained by two studies (Demirtas and
Pechenizkiy, 2013), (Hayran and Sert, 2017). Erkin Demirtas and Mykola
Pechenizkiy (Demirtas and Pechenizkiy, 2013) gathered movie and product
reviews from popular e-commerce and social web sites. The set is gathered from
a popular Turkish cinema web page 777www.beyazperde.com with over 5K positive
and 5K negative comments in the Turkish language. Reviews are marked in scale
from 0 to 5 by the web users. The study considered a review sentiment positive
if the rating is equal to or bigger than 4, and negative if it is less or
equal to 2 as the study does. The products review dataset is taken from an
online e-commerce web page. This benchmark dataset consists of reviews
regarding four product categories. The categories of e-commerce dataset are
dvd, books, electronics and kitchen. Likewise, reviews are marked in the range
from 1 to 5. Each category has equally 700 positive and 700 negative reviews,
where average rating of negative reviews is 2.27 and of positive reviews is
4.5. The dataset has been evaluated in another study (Yildirim, 2020) as well.
Finally we have organized all the datasets for training, development and test
set as described in Table 1

Table 1. Dataset for Sentiment Analysis Sentiment Dataset | # of total comments | Twitter | E-commerce (with Movie) | #Pos | #Neg
---|---|---|---|---|---
Training Set | 40830 | 22000 | 18830 | 20340 | 20490
Development Set | 9000 | 4900 | 4100 | 4445 | 4556
Test Set | 9000 | 4900 | 4100 | 4391 | 4609

The sentiment model was fine-tuned with max-seq-length of 128, learning rate
of 2e-5, number of epoch optimized is 3 888https://huggingface.co/savasy/bert-
base-turkish-sentiment-cased. The performance of the model is reported as
shown in Table 2. We listed other studies performance that used the same
e-commerce and Twitter datasets. We only report their best performance for a
proper comparison. We got 93.12% and 96.69% accuracy for e-commerce and
Twitter dataset respectively. When we merged two datasets for training and
test phase, we got 94.77% accuracy .

Table 2. Sentiment Analysis Performance Evaluation Model | Dataset | Results
---|---|---
Yildirim (2020) | e-commerce | 89.26
Demirtas (2013) | e-commerce | 85.1
Sert (2017) (Hayran and Sert, 2017) | Twitter | 80.5
Our BERT model* | e-commerce | 93.12
Our BERT model* | Twitter | 96.69
Our BERT model* | e-commerce + Twitter | 94.77

###  3.3. Named-Entity Recognition

The fine-tuned model can be found under the link
999https://huggingface.co/savasy/bert-base-turkish-ner-cased Named-Entity
Recognition aims to capture named entities in a text, such as locations,
person. Our NER datasets follow BIO (Beginning, Inside, Outside) data format
which is a common NER tagging format for tokens. The B refers to the first
word of an entity and the I corresponds to the remaining words of the same
entity. Tokens tagged with O means are not part of an entity. B and I tags are
followed by an entity category such as B-PER, I-PER, B-ORG. NER problem is
considered a multi-class token classification task. There are a couple of
named-entity data formats: Enamex, Timex, Numex. Enamex, which is used in our
system, has three entity types: Person, location, organization. NUMEX format
is for numerical entities such as money and percent and TIMEX format is used
for temporal entities such as date and time. Two different datasets annotated
in Enamex format have been used for an evaluation in our fine-tuning phase as
shown in Table 3. The first one is Turkish NER dataset obtained from the study
WikiANN dataset (Pan et al., 2017) which includes 282 different language
dataset compiled from Wikipedia. Second one which is also annotated with POS,
LOC, and ORG is shared by NLP community (Schweter, 2020; Yeniterzi, 2011)

Table 3. NER Dataset Model | Dataset | # of token | # PER | # LOC | # ORG
---|---|---|---|---|---
BERT Model | WikiANN | 340K | 17.2K | 19.1K | 15.8K
BERT Model | Second Dataset | 385K | 16.6K | 10.8K | 10.1K

The initial studies of Turkish NER have started in 1990 (Yeniterzi et al.,
2018). Hidden-markov models, Conditional Random Fields, Bayesian Models and
Neural Networks are among the approaches applied to the NER solution. We
listed NER performances of some Turkish studies (Yeniterzi, 2011; Şeker and
Eryiğit, 2012; Pan et al., 2017; Küçük et al., 2014; Demir and Özgür, 2014;
Şeker and Eryiğit, 2017; Tur et al., 2003) that especially uses Enamex
datasets as shown in Table 4. Bert-based fine-tuning model gained comparable
results of 95.0 F1 and 95.2 F1 score for two datasets as reported in the
table. The study (Pan et al., 2017) still has the highest score of 96.9 F1 by
leveraging cross-lingual language models. They developed a cross-lingual name
tagging and linking framework for 282 different languages that exist in
Wikipedia. The framework could identify name mentions and assign a coarse-
grained type to the mentions. It applies Bi-directional Long Short-Term Memory
and Conditional Random Fields (CRFs) network as the underlying learning model
for the name tagger for each language. It also links the captured names to an
English Knowledge Base.

Table 4. Performance Evaluation of NER Models Model | Dataset | F1 Score |
---|---|---|---
Pan (2017) | WikiANN | 96.9 |
Kucuk (2014) | Twitter data | 72.61 |
Yeniterzi (2018) | Community Set | 88.94 |
Demir (2014) | by (Tur et al., 2003) | 91.85 |
Seker (2012) | by (Tur et al., 2003) | 91.94 |
Seker (2017) | by (Tur et al., 2003) | 92.00 |
Tur (2003) | General News | 92.73 |
Our BERT Model* | WikiANN | 95.0 |
Our BERT Model* | Community Set | 95.2 |

###  3.4. Question Answering

The fine-tuned model can be found under the
link101010https://huggingface.co/savasy/bert-base-turkish-squad

Turkish question answering dataset, namely TQuAD
111111https://github.com/TQuad/turkish-nlp-qa-dataset, has been prepared in
the form of  Stanford Question Answering Dataset (SQuAD) which is a well-known
English QA dataset. The SQuAD includes questions on a set of Wikipedia
articles, where the answer to a question is a span of a text from the
corresponding reading passage. While SQuAD 1.0 format does not cover features
that the question might be unanswerable, this feature was later added on the
SQuAD 2.0 version. The TQuAD dataset currently follows SQuAD 1.0. The examples
have been manually prepared on Turkish & Islamic Science History within the
scope of Teknofest 2018 Artificial Intelligence competition, and it has been
shared with the NLP community. Training and evaluation set have 681 and 72
titles, 2232 and 275 Paragraphs, and 8308 and 892 Q&As respectively. The model
has been fine-tuned with the default BERT model settings such as learning rate
of 0.00003, maximum sequence length of 384. We observed that model
successfully fine-tuned the base pre-trained model with the exact match score
of 62.55 and F1 score of 80.48. For English Language we can see above 95.0 F1
score and we have enough room for enhancing the Question Answering model by
increasing data size and hyper-parameter optimization. Due to the fact that
there is no study using SQuAD-like dataset in Turkish, we are not able to
compare our model to any Turkish study.

###  3.5. Text Classification

The fine-tuned model can be found under the link 121212savasy/bert-turkish-
text-classification. Text classification downstream task has been fine-tuned
with two datasets TTC-4900 (Yildirim, 2014) and TTC-3600 (Kılınç et al.,
2017). TTC-4900 is originally prepared by a Turkish NLP Group
131313http://www.kemik.yildiz.edu.tr and shared with community via repository
141414http://www.kaggle.com/savasy/ttc4900. The data has been already used by
many studies (Yildirim, 2014). The dataset has seven coarse-grained
categories: world, economy, culture, health, politics, sport, and technology.
The dataset equally consists of 700 articles under each category. TTC-3600
dataset has been shared by the study (Kılınç et al., 2017) with a link
151515https://github.com/denopas/TTC-3600. It has six categories; economy,
culture, health, politics, sports, technology. Each category has 600 examples.
The BERT-based model has been fine tuned with these two datasets in the same
settings with the original BERT model and reported as shown in Table 5. As the
table indicates, our fine-tuned models surpassed the traditional approaches
employed in the other studies (Kılınç et al., 2017; Yildirim, 2014) in terms
of F1 score by achieving 93.4 and 92.2.

Table 5. Text Classification Performance Model | Dataset | Size | Performance
---|---|---|---
Yildirim (2014) | TTC 4900 7 Classes | 4900 | 90.0
Kilinc (2017) | TTC 3600 6 Classes | 3600 | 91.3
Our Fine-tuned BERT* | TTC 4900 7 Classes | 4900 | 93.4
Our Fine-tuned BERT* | TTC 3600 6 Classes | 3600 | 92.2

##  4\. Final Remarks

We evaluated four different downstream tasks for the Turkish language. We
gathered several datasets used and shared by the community. Other than
Question Answering dataset, the remaining three datasets have been already
used and evaluated by other studies. We compared our fine-tuned models with
those studies who used the same datasets as shown in Table 2, Table 4, and
Table 5, which suggests that the BERT-based models outperformed many baseline
approaches implemented by other studies.

We observed some well-known facts that when we transfer the model to another
domain, we obtain less performance. For instance, we measure the sentiment
prediction scores ranging from 95.0 F1 to 80.0 F1. But when we use the same
data from the same domain, the performance improves. As we keep fine-tuning
the process fed by that particular domain, the model gets again over 95.0 % F1
score. The Named-Entity Recognition model follows ENAMEX format where three
coarse-grained entity types are considered: PER, LOC, ORG. We want to train
NER models for the problems of NUMEX and TIMEX format. Our preliminary studies
indicated that we can get high performances around 93 F1 with these formats,
which is not reported here. The study (Pan et al., 2017) has an impressive F1
score of 96.9 and better than our models.

Question Answering model cannot cover unanswerable questions, which means no-
answer option is not available and always returns a most likely answer, which
is a standard for SQuAD 1.0. We will extend our model to SQuAD 2.0 by adding a
”no-answer” option. When comparing our QA study to English counterpart, one
can say that both the size of Turkish dataset and the performance of the model
needs to be improved. In the last decades, Text Classification problem has
been successfully solved by traditional models, TF-IDF, or word embeddings
(Yildirim, 2014; Kılınç et al., 2017) in Turkish. We showed that the BERT-
based model has outperformed the traditional approaches for both seven-class
and six-class classification problems. In the future we will fine-tune a model
to cover more fine-grained classes.

There exist many other NLU problems to be addressed in the Turkish Language as
one of the less studied languages. We also plan to leverage other transformer
models such as RoBERTa or ELECTRA to improve the performance of NLU tasks and
solve other types of NLU problems such as anaphora resolution, summarization,
relation extraction in Turkish.

##  5\. Conclusion

The recent studies showed that transformer-based models got state-of-the-art
results in many downstream NLP tasks. In this work, we have leveraged them for
four different downstream NLU tasks: Named-Entity Recognition, Sentiment
Analysis, Question Answering and Text Classification. The experiments showed
that the fine-tuned models significantly outperformed other existing baseline
approaches in the Turkish Language as reported in the study. We publicly
released these four fine-tuned models in reproducibility and in supporting
other researchers. This study is counted to be the first successful attempt to
apply the transformer-based model to the Turkish language by fine-tuning the
BERT-based model. In the future, we want to leverage other types of language
models and cover more challenging down-streams tasks such as coreference
resolution or text summarization for the Turkish language.

## References

  * (1)
  * Bojanowski et al. (2016) Piotr Bojanowski, Edouard Grave, Armand Joulin, and Tomas Mikolov. 2016.  Enriching Word Vectors with Subword Information.  _arXiv preprint arXiv:1607.04606_ (2016).
  * Clark et al. (2020) Kevin Clark, Minh-Thang Luong, Quoc V Le, and Christopher D Manning. 2020.  Electra: Pre-training text encoders as discriminators rather than generators.  _arXiv preprint arXiv:2003.10555_ (2020).
  * Dai et al. (2019) Zihang Dai, Zhilin Yang, Yiming Yang, Jaime G. Carbonell, Quoc V. Le, and Ruslan Salakhutdinov. 2019\.  Transformer-XL: Attentive Language Models Beyond a Fixed-Length Context.  _CoRR_ abs/1901.02860 (2019).  arXiv:1901.02860 http://arxiv.org/abs/1901.02860
  * Demir and Özgür (2014) H. Demir and A. Özgür. 2014.  Improving Named Entity Recognition for Morphologically Rich Languages Using Word Embeddings. In _2014 13th International Conference on Machine Learning and Applications_. 117–122.
  * Demirtas and Pechenizkiy (2013) Erkin Demirtas and Mykola Pechenizkiy. 2013.  Cross-Lingual Polarity Detection with Machine Translation. In _Proceedings of the Second International Workshop on Issues of Sentiment Discovery and Opinion Mining_ (Chicago, Illinois) _(WISDOM ’13)_. Association for Computing Machinery, New York, NY, USA, Article 9, 8 pages.  https://doi.org/10.1145/2502069.2502078
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding. In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_. Association for Computational Linguistics, Minneapolis, Minnesota, 4171–4186.  https://doi.org/10.18653/v1/N19-1423
  * Hayran and Sert (2017) A. Hayran and M. Sert. 2017\.  Analysis on Microblog Data based on Word Embedding and Fusion Techniques. In _2017 25th Signal Processing and Communication Application Conference (SIU)_. tbd.  https://doi.org/tbd
  * Joshi et al. (2019) Mandar Joshi, Omer Levy, Luke Zettlemoyer, and Daniel Weld. 2019\.  BERT for Coreference Resolution: Baselines and Analysis. In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_. Association for Computational Linguistics, Hong Kong, China, 5803–5808.  https://doi.org/10.18653/v1/D19-1588
  * Küçük et al. (2014) Dilek Küçük, Guillaume Jacquet, and Ralf Steinberger. 2014.  Named Entity Recognition on Turkish Tweets. In _Proceedings of the Ninth International Conference on Language Resources and Evaluation (LREC’14)_. European Language Resources Association (ELRA), Reykjavik, Iceland, 450–454.  http://www.lrec-conf.org/proceedings/lrec2014/pdf/380_Paper.pdf
  * Kılınç et al. (2017) Deniz Kılınç, Akın Özçift, Fatma Bozyigit, Pelin Yıldırım, Fatih Yücalar, and Emin Borandag. 2017.  TTC-3600: A new benchmark dataset for Turkish text categorization.  _Journal of Information Science_ 43, 2 (2017), 174–185.  https://doi.org/10.1177/0165551515620551
  * Lample and Conneau (2019) Guillaume Lample and Alexis Conneau. 2019.  Cross-lingual Language Model Pretraining.  _CoRR_ abs/1901.07291 (2019).  arXiv:1901.07291 http://arxiv.org/abs/1901.07291
  * Lan et al. (2020) Zhenzhong Lan, Mingda Chen, Sebastian Goodman, Kevin Gimpel, Piyush Sharma, and Radu Soricut. 2020\.  ALBERT: A Lite BERT for Self-supervised Learning of Language Representations.  _ArXiv_ abs/1909.11942 (2020).
  * 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.  _CoRR_ abs/1907.11692 (2019).  arXiv:1907.11692 http://arxiv.org/abs/1907.11692
  * Mikolov et al. (2013) Tomas Mikolov, Kai Chen, Gregory S. Corrado, and Jeffrey Dean. 2013\.  Efficient Estimation of Word Representations in Vector Space.  _CoRR_ abs/1301.3781 (2013).
  * Pan et al. (2017) Xiaoman Pan, Boliang Zhang, Jonathan May, Joel Nothman, Kevin Knight, and Heng Ji. 2017\.  Cross-lingual Name Tagging and Linking for 282 Languages. In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_. Association for Computational Linguistics, Vancouver, Canada, 1946–1958.  https://doi.org/10.18653/v1/P17-1178
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher Manning. 2014.  Glove: Global Vectors for Word Representation. In _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_. Association for Computational Linguistics, Doha, Qatar, 1532–1543.  https://doi.org/10.3115/v1/D14-1162
  * Peters et al. (2018) Matthew E. Peters, Mark Neumann, Mohit Iyyer, Matt Gardner, Christopher Clark, Kenton Lee, and Luke Zettlemoyer. 2018.  Deep contextualized word representations.  _CoRR_ abs/1802.05365 (2018).  arXiv:1802.05365 http://arxiv.org/abs/1802.05365
  * Radford et al. (2019) Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2019\.  Language Models are Unsupervised Multitask Learners.
  * Raffel et al. (2019) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J. Liu. 2019.  Exploring the Limits of Transfer Learning with a Unified Text-to-Text Transformer.  _ArXiv_ abs/1910.10683 (2019).
  * Schweter (2020) Stefan Schweter. 2020\.  BERTurk - BERT models for Turkish.  https://doi.org/10.5281/zenodo.3770924
  * Şeker and Eryiğit (2012) Gökhan Akın Şeker and Gülşen Eryiğit. 2012\.  Initial explorations on using CRFs for Turkish named entity recognition. In _Proceedings of COLING 2012_. 2459–2474.
  * Tur et al. (2003) Gokhan Tur, Dilek Hakkani-Tur, and Kemal Oflazer. 2003\.  A statistical information extraction system for Turkish.  _Natural Language Engineering_ 9 (06 2003), 181 – 210.  https://doi.org/10.1017/S135132490200284X
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. 2017\.  Attention Is All You Need.  _CoRR_ abs/1706.03762 (2017).  arXiv:1706.03762 http://arxiv.org/abs/1706.03762
  * Wang et al. (2018) Alex Wang, Amanpreet Singh, Julian Michael, Felix Hill, Omer Levy, and Samuel Bowman. 2018\.  GLUE: A Multi-Task Benchmark and Analysis Platform for Natural Language Understanding. In _Proceedings of the 2018 EMNLP Workshop BlackboxNLP: Analyzing and Interpreting Neural Networks for NLP_. Association for Computational Linguistics, Brussels, Belgium, 353–355.  https://doi.org/10.18653/v1/W18-5446
  * Wu et al. (2016) Yonghui Wu, Mike Schuster, Zhifeng Chen, Quoc V. Le, Mohammad Norouzi, Wolfgang Macherey, Maxim Krikun, Yuan Cao, Qin Gao, Klaus Macherey, Jeff Klingner, Apurva Shah, Melvin Johnson, Xiaobing Liu, Lukasz Kaiser, Stephan Gouws, Yoshikiyo Kato, Taku Kudo, Hideto Kazawa, Keith Stevens, George Kurian, Nishant Patil, Wei Wang, Cliff Young, Jason Smith, Jason Riesa, Alex Rudnick, Oriol Vinyals, Greg Corrado, Macduff Hughes, and Jeffrey Dean. 2016.  Google’s Neural Machine Translation System: Bridging the Gap between Human and Machine Translation.  _CoRR_ abs/1609.08144 (2016).  arXiv:1609.08144 http://arxiv.org/abs/1609.08144
  * Yang et al. (2019) Zhilin Yang, Zihang Dai, Yiming Yang, Jaime Carbonell, Ruslan Salakhutdinov, and Quoc Le. 2019\.  XLNet: Generalized Autoregressive Pretraining for Language Understanding.  (06 2019).
  * Yeniterzi (2011) Reyyan Yeniterzi. 2011\.  Exploiting Morphology in Turkish Named Entity Recognition System. In _Proceedings of the ACL 2011 Student Session_. Association for Computational Linguistics, Portland, OR, USA, 105–110.  https://www.aclweb.org/anthology/P11-3019
  * Yeniterzi et al. (2018) Reyyan Yeniterzi, Gökhan Tür, and Kemal Oflazer. 2018\.  _Turkish Named-Entity Recognition_.  Springer International Publishing, Cham, 115–132.  https://doi.org/10.1007/978-3-319-90165-7_6
  * Yildirim (2014) Savaş Yildirim. 2014\.  A Knowledge-Poor Approach to Turkish Text Categorization. In _Proceedings of the 15th International Conference on Computational Linguistics and Intelligent Text Processing - Volume 8404_ (Kathmandu, Nepal) _(CICLing 2014)_. Springer-Verlag, Berlin, Heidelberg, 428–440.  https://doi.org/10.1007/978-3-642-54903-8_36
  * Yildirim (2020) Savaş Yildirim. 2020\.  _Comparing Deep Neural Networks to Traditional Models for Sentiment Analysis in Turkish Language_.  Springer Singapore, Singapore, 311–319.
  * Zhang et al. (2020) Zhuosheng Zhang, Hai Zhao, and Rui Wang. 2020.  Machine Reading Comprehension: The Role of Contextualized Language Models and Beyond.  _ArXiv_ abs/2005.06249 (2020).
  * Şeker and Eryiğit (2017) Gökhan Şeker and Gülşen Eryiğit. 2017.  Extending a CRF-based named entity recognition model for Turkish well formed text and user generated content1.  _Semantic Web_ 8 (01 2017), 1–18.  https://doi.org/10.3233/SW-170253


# A Cacophony of Echoes from daily monitoring of the Crab Pulsar at Jodrell
Bank

Thierry Serafin Nadeau David A. Dunlap Department of Astronomy and
Astrophysics, 50 St. George Street, Toronto, ON M5S 3H4, Canada  Canadian
Institute for Theoretical Astrophysics, 60 St. George Street, Toronto, ON M5S
3H8, Canada  Marten H. van Kerkwijk David A. Dunlap Department of Astronomy
and Astrophysics, 50 St. George Street, Toronto, ON M5S 3H4, Canada  Cees G.
Bassa ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk
4, 7991 PD, Dwingeloo, The Netherlands  Ben W. Stappers Jodrell Bank Centre
for Astrophysics, Department of Physics and Astronomy, University of
Manchester, Manchester M13 9PL, UK  Mitchell B. Mickaliger Jodrell Bank Centre
for Astrophysics, Department of Physics and Astronomy, University of
Manchester, Manchester M13 9PL, UK  Andrew G. Lyne Jodrell Bank Centre for
Astrophysics, Department of Physics and Astronomy, University of Manchester,
Manchester M13 9PL, UK

###### Abstract

Using archival data from the 42 foot telescope at the Jodrell Bank
Observatory, we produce daily stacks of aligned giant pulses for the Crab
pulsar, to study changes to the daily profiles between April 2012 to December
2016. From these, we identify echoes, where intervening material away from the
line of sight causes pulsed emission to be redirected towards the observer,
with delay corresponding to the increased distance of travel, resulting in
additional profile components. These observations show that such echoes may be
far more common than implied by the previous rate of detections. All the
observed echoes are consistent with approaching zero-delay at their closest
approach to the normal giant pulse emission. This indicates that the
structures responsible for producing these events must be highly anisotropic,
with typical lengths greater than $\sim 4\textrm{AU}$, typical widths on the
sky of $\sim 0.1\textrm{AU}$ and typical depths of $\sim 5\textrm{AU}$, given
the previously observed electron densities of the nebular filaments, on the
order of 1000 cm-3. This suggests that these inhomogeneities are likely to be
offshoot substructure from the larger nebular filaments of the Crab nebula.

Pulsars (1306)— Radio Bursts (1339) — Interstellar Scattering (854) — Pulsar
Wind Nebulae (2215)

##  1 Introduction

The Crab pulsar is among the 11 pulsars known to exhibit brief flashes of
extraordinarily bright radio emission known as “giant pulses” (Kuzmin, 2007).

These giant pulses do not occur every rotation and are usually confined a
narrow phase window, often trailing a location where more normal pulsed
emission is seen (e.g., the precursor component in the Crab). For the Crab,
these are the Main Pulse (MP) and Interpulse (IP), which dominate its overall
pulse profile (Hankins & Eilek, 2007). The Crab’s giant pulses are short in
duration, of the order of a few micro-seconds or less, and have been found to
consist of nanosecond-duration “nano shots” (Hankins et al., 2003).

As pulsar radio emission travels to us from its source, it is affected by its
interaction with the ionized plasma of the interstellar medium, or ISM. Hence,
we observe the convolution of the intrinsic pulsar signal with the Impulse
Response Function (IRF) of the intervening material, both via dispersion,
which can be corrected for, and by multipath scattering, which cannot (see
McKee et al. (2018) for a long-term perspective of scattering of the Crab
pulsar). For giant pulses, given their narrow width, the observed signal is
dominated by the IRF, and thus they serve as excellent probes of interstellar
and nebular structures along the line of sight.

As it turns out, on top of being one of the few pulsars with giant pulses, the
Crab is also one of the few pulsars known to sometimes exhibit echoes in its
pulse profile: copies of the profile which are delayed in time relative to the
regular pulse profile. For the Crab, a particularly strong event happened in
1997, which lasted on the order of $\sim\\!100{\rm\;d}$ (Backer et al., 2000;
Lyne et al., 2001). This event made clear that the echoes are caused by the
line of sight passing near structures in the nebula, with the amplitude
increasing and delay decreasing as the line of sight and the structures
approach one another. However, the mechanism underlying the echoes remains
uncertain. Originally, Smith & Lyne (2000) and Backer et al. (2000) proposed
that it was specular reflections off the surface of an intervening ionized
bubble in the Crab nebula. Later, Graham Smith et al. (2011) suggested that
refraction within these clouds was a more likely mechanism.

After the identification of the 1997 event with an echo, similar echoes have
been identified retroactively, in 1974 (Lyne & Thorne, 1975), 1992 and 1994
(Lyne et al., 2001), as well as in observations since (e.g., Driessen et al.
2019). These further echoes are less strong and seem to persist for shorter
periods of time compared to the 1997 event, with timescales of a few days days
to a few weeks (Crossley et al., 2007; Driessen et al., 2019). The delay seen
in the echoes also varies significantly, with Crossley et al. (2007) finding
delays spanning 40 to $100{\rm\,\mu s}$, while Driessen et al. (2019) and
Graham Smith et al. (2011) find delays of 1.5 and $5{\rm\,ms}$, respectively.

As a result of their typically brief duration and relatively irregular
occurrence, echoes are unlikely to be consistently observed in their entirety
with random observations of the Crab pulsar. Instead, regular observations
over longer periods of time are required.

Here, we use Crab observations with the 42 foot telescope at the Jodrell Bank
Observatory for this purpose. We focus on the giant pulses, which, given their
short intrinsic duration and high intensity, serve as very good probes of the
Impulse Response Function of the structures causing the echoes, as long as one
ensures that they are carefully aligned before averaging, so that one avoids
the smearing induced by the fact that giant pulses occur randomly within
$\sim\\!1\%$ of pulse phase.

We describe the observations and the selection of giant pulses in Sect. 2. In
Sect. 3, we describe how they are aligned and averaged into daily stacks, to
gain sufficient signal-to-noise to detect also fainter echoes (with the daily
timescale still well below the timescale on which echoes evolve). In Sect. 4,
we show the results and analyse identified echoes. We discuss the implied
locations and physical properties of the scattering regions in Sect. 5, and
the ramifications of the work in Sect. 6.

##  2 Observations and Giant Pulse Selection

Figure 1:  Identifying genuine giant pulses. Left: Peak signal-to-noise ratio
in the redispersed versus the dedispersed time stream for all triggers in our
data set. The red and green dashed lines indicate the left and upper
boundaries of our selection criterion, removing pulses that are too faint or
that are likely due to broadband RFI – which generally has higher signal-to-
noise in the redispersed timestream. Middle and Right Panels: Dedispersed and
redispersed timestreams for the three triggers marked with blue dots. One sees
that the first two are better aligned in the redispersed timestream,
suggesting they are due to RFI. In contrast, the third trigger is clearly
aligned across frequency, and thus likely is a giant pulse.

Since 1984, the Crab pulsar has been monitored with the 42 foot telescope at
the Jodrell Bank Observatory nearly any time it is visible, mostly to track
the evolution of its spin frequency (Lyne et al., 2015). This work is based on
1721 days of observations, from 2012 April 16 to 2017 January 1 (MJD
56033–57754).

The observations used were obtained with the COBRA2 backend, installed in
2012, which records Nyquist sampled, dual circular polarization, complex
voltages for 10 MHz of bandwidth centered at 610 MHz, with a strong dropoff in
signal outside the central 5 MHz due to an RF filter. The resulting complex
voltages were coherently dedispersed at the nominal dispersion measure (DM) of
the Crab pulsar of $\sim\\!56.8{\rm\,pc/cm^{3}}$. The data were channelized
into 0.5 MHz channels using dspsr (van Straten & Bailes, 2011) and folded
using the regularly updated Jodrell Bank Crab ephemerides (Lyne et al., 1993).
Pulse profiles for individual rotations were processed by masking the outer 4
frequency channels and averaging each profile in frequency and polarization to
512 pulse phase bins. The standard psrchive (Hotan et al., 2004) metrics were
used to define the off-pulse region, and its mean and standard deviation used
were to select those pulse rotations where the maximum flux value of the
entire profile exceeded the noise by $6\sigma$. For those profiles, pulse
profiles with full frequency and polarization information were stored at 8192
pulse phase bins, where each bin contains 2 or 3 channelized time samples.

Since the selection is based on a simple signal-to-noise (SN) threshold, the
archived profiles contain many false positives due to radio frequency
interference (RFI). One common type of RFI is narrow-band, being more or less
continuously present in the top two frequency channels (614–615 MHz). Since
this RFI dominates those channels causing both false positives and substantial
extra noise in genuine pulses, we mask these channels in our analysis.

A second common form of RFI is impulsive and broadband (for an example, see
Figure 1). In the dedispersed time stream, it results in profiles that last
roughly 0.01s, sweeping up in frequency, i.e, with the time delay expected
from dedispersion over the 5-MHz frequency band for the pulsar (DM) of
$\sim\\!56.8{\rm\,pc/cm^{3}}$.

Since genuine giant pulses are narrow in the dedispersed time stream, we can
distinguish those from the impulsive RFI by comparing the SN as measured in
the dedispersed time stream with the SN inferred after redispersing the data
to its original state: genuine giant pulses and RFI will have higher SN in the
dedispersed and redispersed time streams, respectively. In Figure 1, we show
the SN with and without redispersion for all triggers, along with a few
examples. The real giant pulses are confined to a band with a range in
dedispersed peak SN but low redispersed SN.

For our selection of giant pulses, we start by removing all low-SN triggers,
with a dedispersed peak SN of less than $5\sigma$ (thus removing pulses
triggered or greatly affected by the narrow-band RFI). Next, we exclude all
triggers above the green cutoff shown in Figure 1, which is flat (redispersed
S/N of $4.5\sigma$) before transitioning to a slope of 1 for the brightest
pulses (where the numbers were chosen empirically to exclude as many false
triggers as possible while being sure to include pulses that visual inspection
showed were clearly genuine). In the end, these cuts removed 2256672 of
3038254 triggers, or about 74.3%.

Even with this initial SN based selection, the distribution over pulse phase
shows that there are still many profiles incorrectly identified as as Crab
giant pulses (see Figure 2). In order to reduce the false-positive rate
further, we ignore all triggers that are more than 1.5% in phase away from the
Main pulse and Interpulse phase windows. For this purpose, we first convert
times to pulsar phases using `tempo2` (Hobbs et al., 2006) along with the
Jodrell Bank Crab ephemerides (Lyne et al., 1993) for each sidereal day. We
find that in the resulting phase distributions, the main and interpulse wander
slowly with time as the ephemeris loses its predictive power, and shows jumps
when the timing ephemeris is updated. For our purposes, though, this is not
very important: for almost all days, we can use the distribution of triggers
in phase to reliably determine the offsets of the Main pulse and Interpulse
windows, and only a few real pulses end up being removed by our phase gating.
Lastly, during phase gating, we mask out times in which there is an excess of
triggers across all phases, implying a large rate of false positives also in
the pulse gates. After the phase selection, we are left with a total of 733663
triggers.

We can estimate the remaining false-positive rate from the number of profiles
which are out of phase with the giant pulses after our SN selection but not at
bad times (437954 triggers). Given that these cover 94% in phase, we infer
that in each of the Main and Interpulse phase windows there will be about
14000 false positives, or a false-positive rate of 3.8% (2.2% and 16% of Main
Pulse and Interpulse, respectively). This of course varies from day to day
with the 1$\sigma$ range around the median being $3.1^{+10}_{-0.7}\%$.

Figure 2:  Phase distribution of all triggers (top), after removing low
signal-to-noise ones (middle), and after clipping against broadband RFI
(bottom), using constraints shown in Figure 1. One sees that a large majority
of archival triggers profiles are not genuine giant pulses.

##  3 Pulse Profiles

Figure 3:  Pulse stacks, before (top) and after (middle) alignment. Stacks are
shown for three days, with low, moderate and high number of giant pulse
detections, with the resulting weighted profiles below. For all days,
correcting for jitter in the exact phase at which a pulse occurs leads to
clear improvements, with the most drastic effect seen in days with low
scattering, when the phase jitter is large compared to the scattering time.
Fits to the aligned and averaged stack profiles are also shown (bottom) with
the residuals from the fits below, for both the exponential tail (green dot-
dashed line) and the modified exponential tail (blue dashed line) functions.

For different days, the removal of false triggers leaves us with greatly
varying numbers of pulses, ranging from tens to thousands. Generally, fewer
pulses are detected on days that the pulses have very long scattering tails as
their peak fluxes thus less frequently reach the trigger threshold. On days
when the Sun is close to the Crab on the sky, the number of detections is also
reduced due to the decreased sensitivity of the telescope. More accurate pulse
profiles are obtained by averaging pulses, but before doing that, we must take
into account that pulses do not arrive at exactly the same phase, but rather
randomly occur in their phase windows, i.e., up to a few hundreds of $\mu s$
apart (Figure 3), a phenomenon intrinsic to the emission known as jitter.

In order to determine the individual pulse phases, we fit each profile with a
model consisting of a (narrow) Gaussian convolved with an exponential
scattering tail. We then use the positions inferred from these fits to shift
all pulses to a common phase (multiplying with a phase ramp in the Fourier
domain to most easily allow for sub-pixel shifts). Next, we create summed
pulse profiles for each sidereal day, optimally weighting them by normalizing
each pulse by the square of the off-pulse rms noise.

It was found, however, that the fitting process is not entirely reliable for a
small portion of the individual profiles, with the fits either failing
entirely or resulting in visibly incorrect pulse positions. Furthermore,
fitting the same exponential tail model to the resulting average profiles, it
becomes clear that this is not always a good model: the scattering tail is
often sharper than an exponential tail allows for, falling more quickly near
the start and more slowly in the tail. Better results were obtained with a
model consisting of a Gaussian convolved with a “modified exponential tail”,
of the form $\exp((-t/\tau)^{\gamma})$, where $\gamma$ determines the
sharpness. With this added parameter, the observed scattering tails are much
more consistently reproduced (see Fig. 3). Values of $\gamma$ between 1 and
0.25 were obtained, with lower scattering times typically corresponding to
lower $\gamma$, and larger scattering times corresponding to larger $\gamma$.

To align the profiles more reliably, we iterate: we fit the first average we
constructed with the improved modified exponential tail model, and then use
those fits as a matched filter for the individual pulse profiles of the
corresponding sidereal day. We then obtain improved positions by fitting the
resulting peaks with a Gaussian, use those positions to align the profiles,
and sum those together as before to construct the final average profiles.

##  4 A Cacophony of Echoes

Fig. Set4. Evolution of the Scattering Profile

Figure 4:  Evolution of the scattering profile. Top: Stack of daily average
giant pulse profiles for the full data set. Middle: Expanded view of the
average profiles indicated by the blue box (MJD 56632 – 56831). Bottom:
Residuals from fits to these average profiles with a narrow Gaussian convolved
with a modified exponential tail (shown on a non-linear colour scale to help
highlight echoes). In the expanded views, the times of closest approach of
echoes are marked by green ticks. For the more prominent echoes, in the
residuals panel parabolic arcs are overlaid, with dashed line style for the
side where the echo was visible, and dotted style for the continuation on the
other side. The complete figure set (12 figures) is available in the online
journal.

The daily average profiles show the strong scattering variations expected from
the Crab’s scattering environment, as is immediately obvious in Figure 4,
where they are combined into a stack. Also visible in the stack, especially
when zoomed in to a finer timescale, are prominent echoes, some with incoming
and some with outgoing arcs.

Inspecting the data set, a number of things stand out. First, strong echoes
often appear clustered together, suggesting a possible causal connection,
perhaps to some nebular structure.

Second, the clearest echoes appear during days with lower scattering. This
might be a signal-to-noise issue: fewer giant pulses are detected on days with
higher scattering, and echoes are more difficult to see against a longer
scattering tail. Given the observed clustering of prominent echoes, however,
it may also be that observed periods of high scattering are at least partly
the result of the presence of many echoes, which blend together to give the
appearance of a longer scattering tail.

Third, during periods of low scattering, many less prominent echoes appear to
be present. Unlike the clearer cases discussed above, these are difficult to
see directly in the stacks but seem convincingly identified as a repeating
pattern of arclet structures in the stacks of the residuals between the daily
average profiles and the fitted scattering tail model (see Fig. 4).

Fourth, the observed echoes trace arcs that are consistent with approaching
zero delay, implying that the structures responsible closely approach or
directly cross the line of sight to the Crab. This is unexpected if the
structures responsible are some kind of “blobs,” as for any roughly spherical,
localized structures directly crossing the line of sight should be the
exception rather than the rule, i.e., the typical observed minimum delay of
the echoes should be nonzero. Instead, it suggests structures that are highly
anisotropic (as projected on the sky), causing an echo that moves along their
long axis, such that as the pulsar moves, the echo will be on the line of
sight when it crosses the structure.

Fifth, the arcs look parabolic, with the delay depending quadratically on time
relative to the zero-delay crossing. This is expected generically for linear
motion relative to a fixed scattering structure, with the curvature depending
on a combination of distance, orientation, and relative velocity (see Sect.
6). We matched parabolae to the most prominent echoes in our data set (see
Fig. 4) and found that the curvatures are not the same, but vary between $6$
and $17{\rm\,\mu s/day^{2}}$, implying differences in distances, orientations
or velocities. We also find that while echoes clustered together in time have
similar curvatures, within $2{\rm\,\mu s/day^{2}}$, they are not identical.
Thus, if these are caused by the same larger-scale structure and thus have the
same distance and velocity, their relative orientations must differ.

Sixth, even the most prominent echoes only rarely, if ever, seem to continue
across zero (where they have parts on both sides with the same curvature).
This suggests that a given structure can typically bend radiation only in one
direction.

Overall, the results show that echoes are more common than previously
realized, and will be apparent when care is taken to correct for jitter in
individual pulses before stacking. They also suggests that a “cacophony of
echoes” is always present and that echoes may at least be partially
responsible for the extent of scattering tails. Furthermore, the structures
responsible are likely highly elongated (as projected on the sky) and can bend
light in only one direction.

##  5 Properties of the Echo Structures

The scattering is thought to occur in the Crab nebula, and is likely related
to the filaments that are seen in optical emission lines. Given that those
filaments are highly elongated, it is not unexpected that structures
associated with them on smaller scales, including those that could cause the
echoes we see, are elongated as well. Below, we check whether our observations
are consistent with this picture, first discussing constraints on the
locations of the scattering structures, and then turning to their shapes.

###  5.1 Locations

An elongated structure in the Crab nebula (i.e., much closer to the pulsar
than to us), which is crossed by the pulsar’s projected trajectory, will
produce an echo with a geometric delay $\tau$ that varies quadratically with
time as,

$\tau(t)=\delta(t)^{2}\frac{d_{ps}}{2c}=\frac{\left(v_{\rm
eff}\cos\psi\right)^{2}}{2cd_{ps}}t^{2}\equiv\eta t^{2},$ (1)

where $\delta(t)$ is the angular offset as a function of time $t$ (relative to
the crossing time), $d_{ps}$ the distance between the pulsar and the screen,
and $\psi$ is the angle between the position angle of the normal to the
structure and the effective velocity $v_{\rm eff}$, and where in the final
part we implicitly defined the curvature $\eta$. The effective velocity is
given by the difference between the pulsar and structure velocities projected
onto the sky,

$\vec{v}_{\rm eff}=\vec{v}_{p,\rm sky}-\vec{v}_{s,\rm sky}.$ (2)

Note that these are just the normal equations for scattering of a thin screen
(e.g., Cordes & Rickett 1998), but with the simplification that the screen is
close to the pulsar.

Figure 5:  Distribution of the measured curvatures for the most prominent
echoes. Overlaid are curves expected to be generated by structures that are
located anywhere in the nebula with random orientations ($d_{ps}\sim{\cal
U}(0.5,2.0){\rm\,pc}$; dotted line), at the edge of the wind nebula with
random orientations ($d_{ps}=0.5{\rm\,pc}$; dashed line), and near the edge of
the nebula with a small range of orientations ($d_{ps}\sim{\cal
N}(0.5,0.1){\rm\,pc}$, $\psi\sim{\cal N}(40.7,0.6){\rm\,deg}$; full line).

The distribution of the curvatures $\eta$ for the observed echoes is shown in
Fig. 5. Most curvatures lie in the range between $\eta=$6 and
$\eta=10{\rm\,\mu s/day^{2}}$, while a few have higher values, yielding an
average curvature of $\bar{\eta}=9.7{\rm\,\mu s/day^{2}}$.

For a given echo arc, the curvature $\eta$ constrains a combination of the
corresponding structure’s distance from the pulsar, its effective velocity,
and its orientation, with the highest curvatures providing the most stringent
constraints. From observations of nebular emission, Martin et al. (2021) found
that the optical filaments roughly extend from 0.5 to $2.0{\rm\,pc}$ (where
the inner boundary corresponds to the edge of the pulsar wind nebula; H. E. S.
S. Collaboration 2020). If the material responsible for producing the fastest
echoes, with $\eta=17{\rm\,\mu s/day^{2}}$ is located at least $0.5{\rm\,pc}$
from the pulsar, it requires that $v_{\rm eff}\gtrsim 145\,\textrm{km/s}$.
This is slightly higher than the pulsar’s velocity of $v_{p}=120\pm
23{\rm\,km/s}$ inferred from astrometry (Kaplan et al., 2008), though still
consistent given the relatively large uncertainties on the pulsar velocity
relative to the center of the nebula and the possible contribution to the
effective velocity of the motion of nebular material.

Scaling $v_{\rm eff}$ to this velocity, an estimate of the distance to any
echo is then,

$d_{ps}=0.85{\rm\;pc}\left(\frac{\eta}{10{\rm\,\mu
s/day^{2}}}\right)^{-1}\left(\frac{v_{\rm
eff}\cos\psi}{145{\rm\,km/s}}\right)^{2},$ (3)

which is an upper limit, since $\cos\psi\leq 1$. From our observed curvature
range of $6\leq\eta\leq 17{\rm\,\mu s/day^{2}}$, we therefore infer a range in
maximum distances of 0.5 to $1.4{\rm\,pc}$. Here, the lower bound matches the
Martin et al. (2021) range by construction, but the fact that the upper
boundary is also within the expected range from the nebular emission gives
some confidence in the estimates.

The distribution in curvature gives further information, as can be seen in
Fig. 5, where we compare it with expectations for three distinct cases. The
first two assume random screen orientations, with structures either placed
uniformly between 0.5 and $1.5{\rm\,pc}$ or at a fixed distance of
$0.5{\rm\,pc}$. Neither matches, as they predict a peak at very low curvature
that is not seen (corresponding to $\cos\psi\simeq 0$, i.e., structures nearly
aligned with the pulsar motion so that the separation between the structure
and the line of sight varies slowly).

In the third case, we sample both screen distance and orientation from narrow
Gaussian distributions. This case reproduces the observed distribution of echo
curvatures much better. As expected from Eq. 3, there is a degeneracy between
the mean distance and mean orientation: given a minimum mean screen distance
$\bar{d}_{ps}>0.5{\rm\,pc}$ and assuming $v_{\rm eff}=145{\rm\,km/s}$, we
constrain the mean orientation to $\bar{\psi}\leq 40.9\arcdeg$; conversely,
given that $\bar{\psi}\geq 0\arcdeg$, we constrain the mean distance to
$\bar{d}_{ps}\leq 0.87{\rm\,pc}$.

Regardless of the precise distance or orientation, the range in both $d_{ps}$
and $\psi$ has to be relatively small. We return to possible interpretations
of this in Sect. 6.

Figure 6: Geometry of an echo, with an on-sky view on the left, and one along
the line of sight on the right. The grey line is the structure responsible for
the echoes (which can only bend light in one direction), and the blue arrow
marks the pulsar’s path. The yellow lines trace the rays of the lensed image
at three different times (with corresponding line-of-sight paths shown in
green on the right). In the left-hand panel, four regions are defined: (A)
Echoes will be seen starting at zero delay until the maximum distance for
which the structure can still bend light towards the line of sight; (B) Echos
will be seen starting at non-zero delay; (C) The pulsar passes at a distance
too far from the structure to see echoes; and (D) The pulsar passes on the
wrong side of the structure to see echoes. Note that the right-hand panel is
strongly compressed along the line of sight; we infer from our observations
that the width (vertical extent here) is far smaller than the depth
(horizontal extent). The orientation angle of the screen relative to the
pulsar motion is $\psi$, $\delta$ is the bending angle, and $\theta$ and
$\beta$ are the angular offsets between the structure and the echo and line of
sight images respectively.  Figure 7: Sketch of the large-scale geometry of
the structures causing echoes, with an on-sky view on the left, and one along
the line of sight on the right. Like in Figure 6, the blue arrow is the path
of the pulsar. Filaments are drawn as thick grey lines, with regions that may
exhibit echoes highlighted in yellow. Short “hairs” along the filaments
represent the ripples or smaller-scale Rayleigh-Taylor instabilities which we
suggest may be directly responsible for the echoes.

###  5.2 Shapes

The strengths and durations of echoes are influenced by the typical electron
densities $n_{e}$ of the structures, as well as their size and shape. Hence,
we can use the echoes and estimates of $n_{e}$ from optical emission lines to
constrain the typical length $L$, width $W$ and depth $Z$ of the lens
structures.

The first qualitative constraint we have is that the structures must be highly
elongated, i.e., $L\gg W$, in order for the echoes to last long and have a
high likelihood of approaching a minimum delay of zero (see Sect. 4). As the
lensed image moves it will be on the part of the lens closest to the pulsar,
moving along with it (see Fig. 6). Thus for any observed echo, one can
estimate a minimum size $\ell_{i}$ that is needed to ensure the lens is
visible for a given duration $\Delta t_{i}$:

$\displaystyle\ell_{i}$ $\displaystyle=$ $\displaystyle v_{\rm eff}\Delta
t_{i}\sin|\psi_{i}|$ (4) $\displaystyle=$ $\displaystyle
0.8{\rm\,AU}\,\sin|\psi_{i}|\left(\frac{\Delta
t_{i}}{10\textrm{days}}\right)\left(\frac{v_{\rm
eff}}{145{\rm\,km/s}}\right).$

Furthermore, for a given structure with total length $L_{i}$, the probability
that the resulting echo will approach a minimum delay $\tau_{\rm
min}=0{\rm\,\mu s}$ is the same as the probability that the line of sight will
intersect the lens as the pulsar moves on the sky (see Fig. 6), giving
$P_{{\rm LOS},i}=L_{i}/(L_{i}+\ell_{i})$. Under the simplifying assumption
that these structures share a typical length $L_{i}\sim L$ and that
$L\gg\ell_{i}$, the probability that out of N observed echoes, all have
$\tau_{\rm min}=0{\rm\,\mu s}$ is then,

$P_{LOS}(N)=\prod_{i=1}^{N}\frac{L}{L+\ell_{i}}\simeq\left(\frac{L}{L+\langle\ell\rangle}\right)^{N}$
(5)

where $\langle\ell\rangle$ is some suitable average of the estimates for
individual echoes.

In order to ensure that such an occurrence is not too rare, we require
$P_{LOS}(N)\geq P_{c}$, where $P_{c}$ is the desired (small) probability that
seeing only zero-crossing echoes was a fluke. Rearranging, this corresponds to
a lower limit to the length,

$L\geq\frac{P_{c}^{1/N}}{1-P_{c}^{1/N}}\langle\ell\rangle.$ (6)

Given that in our $N=22$ observed echoes with measured curvatures, none have
non-zero minimum delays, we infer that with 90% confidence (i.e.,
$P_{c}=0.1$),

$L\geq 4{\rm\,AU}\,\left(\frac{\sin(\bar{\psi})}{0.5}\right)\left(\frac{\Delta
t}{10\textrm{days}}\right)\left(\frac{v_{\textrm{eff}}}{145{\rm\,km/s}}\right).$
(7)

To constrain the width W and depth Z of the structures, we need to make an
assumption about how they create echoes. For the width, we will assume that
they act like lenses, and for the depth that light is bent by gradients in
electron column density.

By conservation of surface brightness, generally the expected magnification
for a lens is $\mu={\rm d}\theta/{\rm d}\beta$, where $\theta$ and $\beta$ are
the echo and source positions relative to the structure, respectively. The
position $\theta$ is constrained to be within the lens, being offset by a
small amount when the echo is bright and approaching the angular half-width
$\omega\equiv W/2d_{ps}$ as it becomes faint. Hence, for a faint image, one
has $\mu\simeq\omega/\beta$. Furthermore, for a faint image, $\omega\ll\beta$
and hence the angular offset between the source and the echo image
$\delta=\beta-\theta\simeq\beta$.

Given the above, we can estimate the physical width of the lens with,

$W\simeq 2\mu\delta d_{ps}\simeq
0.1{\rm\;AU}\left(\frac{\mu}{0.05}\right)\left(\frac{\delta}{1\arcsec}\right)\left(\frac{d_{ps}}{1{\rm\;pc}}\right),$
(8)

where we scaled $\delta$ to a typical value (from Eq. 1 using $\tau\simeq
1{\rm\,ms}$ and $d_{ps}\simeq 1{\rm\,pc}$).

To constrain the depth $Z$, we use that for a given bending angle – in our
case the same as the angular offset $\delta$ implied by an echo’s delay $\tau$
– the required gradient in electron column density $N_{e}$ is,

$\nabla_{x}N_{e}=\frac{2\pi\delta}{\lambda^{2}r_{e}}\simeq 4.5\times
10^{4}{\rm\;cm^{-2}\,cm^{-1}}\left(\frac{\delta}{1\arcsec}\right),$ (9)

where $x$ is the direction along which images form, perpendicular to the
structure, $\lambda\simeq 0.49{\rm\,m}$ is the observing wavelength, and
$r_{e}$ is the classical electron radius.

For a given change in electron density $\Delta n_{e}$, the column density
gradient can be approximated as $\nabla_{x}N_{e}\simeq\Delta n_{e}(Z/W)$, and
hence one can estimate $Z=W\nabla_{x}N_{e}/\Delta n_{e}$. Scaling to the
electron density typical of observed nebular filaments, $\Delta n_{e}\simeq
1000\textrm{cm}^{-3}$ (Osterbrock, 1957), we infer,

$Z=45\,W\left(\frac{\delta}{1\arcsec}\right)\left(\frac{\Delta
n_{e}}{1000\textrm{cm}^{-3}}\right).$ (10)

Thus, we find that the structures responsible for the echoes likely have
sheet-like shapes, with thicknesses $W\approx 0.1{\rm\,AU}$, and sizes
$L\gtrsim 4{\rm\,AU}$ and $Z\approx 5{\rm\,AU}$.

##  6 Ramifications

We have shown that alignment of pulses before stacking allows for very
sensitive observations of echoes. In so doing, we have revealed that echoes
are extremely common.

We find that the echoes are likely caused by sheet-like structures being
viewed edge on, with structures elongated in only one dimension clearly
excluded. The length scales of the sheets, of $\gtrsim\\!4{\rm\,AU}$ on the
sky, are below the resolution of present observations of nebular filaments
($\sim\\!100{\rm\,AU}$ with _HST_), but given that these are thought to be the
result of Rayleigh-Taylor instabilities (Hester, 2008), might well be
substructures of those filaments. Indeed, that the structures are part of a
larger filament is also suggested by the clustering of the curvatures of many
of the echoes, which we found requires sharing similar orientation and
distance.

A possible interpretation is that the structures causing the echoes are nearly
parallel ripples or smaller Rayleigh-Taylor fingers of a larger structure,
perhaps a nebular filament, that happens to be crossing our line of sight.
Physically, the echoes might arise in ionized interfaces of the neutral
interiors of these structures with the hot pulsar wind material. This would
naturally be sheet-like and, when viewed from the side, be able to to bend
light sufficiently to cause echoes. We would also expect these to therefore be
coincident with changes in the dispersion measure, as was observed in the 1997
event reported by Backer et al. (2000) and Lyne et al. (2001).

In this picture, one would expect that the typical curvature of the echoes
would change only gradually, until, as the pulsar moves, it leaves larger
filaments behind and make its way to new ones, with new typical distances and
orientations, and thus different typical curvatures (see Fig. 7). Indeed, from
previous work, we know that the echo curvatures can vary greatly. The 1997
echo evolved with a curvature of $2{\rm\,\mu s/day^{2}}$, much lower than the
curvature of any echo that we see (and hence lasted much longer, though it was
also far more prominent). This is consistent with the Crab now having moved
well past any structure associated to the 1997 event.

Over the time frame observed in this data, we see the 2012 December echo
previously identified in Driessen et al. (2019), but not that seen in 2013
January. For the 2012 December event, our observed delay at $610\text{MHz}$
appear shorter, $\sim\\!700{\rm\,\mu s}$, than that seen at $350{\rm\,MHz}$,
$\sim\\!1350{\rm\,\mu s}$. Backer et al. (2000) showed such differences in
delay may be expected for both refraction and dispersion-based scenarios;
better frequency coverage would be needed to distinguish them. It is not clear
why we do not observe the 2013 January echo seen by Driessen et al. (2019);
perhaps at our higher frequency it could no longer bend light sufficiently.
Conversely, that Driessen et al. (2019) did not see many of our smaller echoes
is not surprising: they looked at folded pulse profiles rather than stacked,
aligned giant pulses and hence would not be able to see our typically fainter
echoes at short delays.

While our stacking of giant pulses has greatly improved the sensitivity to
echoes, significant improvement would come by observing with larger bandwidth.
This would allow one to analyze the frequency dependency of echoes directly,
and also increase the sensitivity to fainter events, thus allowing stacks of
shorter duration and reduced smearing of echo profiles. We plan to pursue this
with daily monitoring of the Crab with the Canadian Hydrogen Intensity Mapping
Experiment (CHIME; CHIME Collaboration et al., 2022).

Acknowledgements: We thank the Toronto scintillometry group for useful
discussions. M.Hv.K. is supported by the Natural Sciences and Engineering
Research Council of Canada (NSERC) via discovery and accelerator grants, and
by a Killam Fellowship. Pulsar research at Jodrell Bank Observatory and the
use of the 42ft telescope is supported by a consolidated grant from the UK
Science and Technology Facilities Council (STFC). Computations were performed
on the Sunnyvale computer at the Canadian Institute for Theoretical
Astrophysics (CITA).

Software: Software: astropy (Astropy Collaboration et al., 2013, 2018, 2022),
numpy (Harris et al., 2020), matplotlib (Hunter, 2007), scipy (Virtanen et
al., 2020), tempo2 (Hobbs & Edwards, 2012).

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Backer et al. (2000) Backer, D. C., Wong, T., & Valanju, J. 2000, ApJ, 543, 740, doi: 10.1086/317150
  * CHIME Collaboration et al. (2022) CHIME Collaboration, Amiri, M., Bandura, K., et al. 2022, ApJS, 261, 29, doi: 10.3847/1538-4365/ac6fd9
  * Cordes & Rickett (1998) Cordes, J. M., & Rickett, B. J. 1998, ApJ, 507, 846, doi: 10.1086/306358
  * Crossley et al. (2007) Crossley, J. H., Eilek, J. A., & Hankins, T. H. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 365, SINS - Small Ionized and Neutral Structures in the Diffuse Interstellar Medium, ed. M. Haverkorn & W. M. Goss, 271.  https://arxiv.org/abs/astro-ph/0612109
  * Driessen et al. (2019) Driessen, L. N., Janssen, G. H., Bassa, C. G., Stappers, B. W., & Stinebring, D. R. 2019, MNRAS, 483, 1224, doi: 10.1093/mnras/sty3192
  * Graham Smith et al. (2011) Graham Smith, F., Lyne, A. G., & Jordan, C. 2011, MNRAS, 410, 499, doi: 10.1111/j.1365-2966.2010.17459.x
  * H. E. S. S. Collaboration (2020) H. E. S. S. Collaboration. 2020, Nature Astronomy, 4, 167, doi: 10.1038/s41550-019-0910-0
  * Hankins & Eilek (2007) Hankins, T. H., & Eilek, J. A. 2007, ApJ, 670, 693, doi: 10.1086/522362
  * Hankins et al. (2003) Hankins, T. H., Kern, J. S., Weatherall, J. C., & Eilek, J. A. 2003, Nature, 422, 141, doi: 10.1038/nature01477
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hester (2008) Hester, J. J. 2008, ARA&A, 46, 127, doi: 10.1146/annurev.astro.45.051806.110608
  * Hobbs & Edwards (2012) Hobbs, G., & Edwards, R. 2012, Tempo2: Pulsar Timing Package, Astrophysics Source Code Library, record ascl:1210.015.  http://ascl.net/1210.015
  * Hobbs et al. (2006) Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655, doi: 10.1111/j.1365-2966.2006.10302.x
  * Hotan et al. (2004) Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302, doi: 10.1071/AS04022
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Kaplan et al. (2008) Kaplan, D. L., Chatterjee, S., Gaensler, B. M., & Anderson, J. 2008, ApJ, 677, 1201, doi: 10.1086/529026
  * Kuzmin (2007) Kuzmin, A. D. 2007, Ap&SS, 308, 563, doi: 10.1007/s10509-007-9347-5
  * Lyne et al. (2015) Lyne, A. G., Jordan, C. A., Graham-Smith, F., et al. 2015, MNRAS, 446, 857, doi: 10.1093/mnras/stu2118
  * Lyne et al. (1993) Lyne, A. G., Pritchard, R. S., & Graham Smith, F. 1993, MNRAS, 265, 1003, doi: 10.1093/mnras/265.4.1003
  * Lyne et al. (2001) Lyne, A. G., Pritchard, R. S., & Graham-Smith, F. 2001, MNRAS, 321, 67, doi: 10.1046/j.1365-8711.2001.03998.x
  * Lyne & Thorne (1975) Lyne, A. G., & Thorne, D. J. 1975, MNRAS, 172, 97, doi: 10.1093/mnras/172.1.97
  * Martin et al. (2021) Martin, T., Milisavljevic, D., & Drissen, L. 2021, MNRAS, 502, 1864, doi: 10.1093/mnras/staa4046
  * McKee et al. (2018) McKee, J. W., Lyne, A. G., Stappers, B. W., Bassa, C. G., & Jordan, C. A. 2018, MNRAS, 479, 4216, doi: 10.1093/mnras/sty1727
  * Osterbrock (1957) Osterbrock, D. E. 1957, PASP, 69, 227, doi: 10.1086/127053
  * Smith & Lyne (2000) Smith, F. G., & Lyne, A. G. 2000, in Astronomical Society of the Pacific Conference Series, Vol. 202, IAU Colloq. 177: Pulsar Astronomy - 2000 and Beyond, ed. M. Kramer, N. Wex, & R. Wielebinski, 499
  * van Straten & Bailes (2011) van Straten, W., & Bailes, M. 2011, PASA, 28, 1, doi: 10.1071/AS10021
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2


11institutetext: Department of Clinical Sciences, Lund University, Sweden,
11email<EMAIL_ADDRESS>
22institutetext: Cognitive Systems, Department of Applied Mathematics and
Computer Science,
DTU, Denmark, 22email<EMAIL_ADDRESS>
33institutetext: Department of Statistics, Lund University, Sweden,
33email<EMAIL_ADDRESS>
# EMPIRICALLY DRIVEN SPLINE BASES FOR FUNCTIONAL PRINCIPAL COMPONENT ANALYSIS
IN 2D

Rani Basna∗ 11 Hiba Nassar ∗ 22 Krzysztof Podgórski∗ 33

###### Abstract

Functional data analysis is typically performed in two steps: first,
functionally representing discrete observations, and then applying functional
methods to the so-represented data. The initial choice of a functional
representation may have a significant impact on the second phase of the
analysis, as shown in recent research, where data-driven spline bases
outperformed the predefined rigid choice of functional representation. The
method chooses an initial functional basis by an efficient placement of the
knots using a simple machine-learning algorithm. The approach does not apply
directly when the data are defined on domains of a higher dimension than one
such as, for example, images. The reason is that in higher dimensions the
convenient and numerically efficient spline bases are obtained as tensor bases
from 1D spline bases that require knots that are located on a lattice. This
does not allow for a flexible knot placement that was fundamental for the 1D
approach. The goal of this research is to propose two modified approaches that
circumvent the problem by coding the irregular knot selection into their
densities and utilizing these densities through the topology of the spaces of
splines. This allows for regular grids for the knots and thus facilitates
using the spline tensor bases. It is tested on 1D data showing that its
performance is comparable to or better than the previous methods.

###### keywords:

Splinets, tensor spline bases, orthonormal bases, binary regression trees

**footnotetext: All authors contributed equally.

##  1 Introduction

The data-driven knot selection algorithm (DDK), see [1], identifies regions in
the common domain of functional data that feature high curvatures and other
significant systematic variability. By strategically placing knots in these
regions, the algorithm effectively captures the systematic variation and is
available as a package on the GitHub repository [2]. In the one-dimensional
domain case, the placement of knots can be directly utilized in the choice of
the splines spanned over selected knots and it has been proved that this
improves the efficiency of the functional representation of the data and the
ensuing statistical analysis. The reason behind the improvement is dimension
reduction in the initial representation of the data achieved by a ‘clever’
choice of the knots. Figure 1 (top) presents the knot selection for the wine
spectra (left) and their derivatives (right), see [3]. Once the knots are
provided, one can build the splinet i.e. a convenient orthonormal spline basis
spanned over these knots, [4]. How the choice of the knots can improve the
analysis of the data is shown in Figure 1 (bottom), where we see the
projection of the data to the first four principal components for equally
spaced knots (left) and for the optimized data-driven knot selection (right).

Figure 1: 1D data and the knot selection. (Top): 1D case, the wine data, the
knots are marked by vertical dashed lines. (Bottom): Projection to the first
four principal component for equally spaced (left) and data-driven, optimally
spaced knots (right).

In 2D, one can utilize tensor product spline bases only, since the knots
reside on the regular lattice knots. On the other hand, if one chooses ‘knots’
to optimize the mean square error (MSE), the selected ‘knots’ are irregular
and they cannot be used directly to build convenient spline bases. In this
note, we discuss how an irregular ‘knot’ selection can be utilized to modify
the approach that allows the use of the tensor bases but at the same time
accounts for the data-driven ‘knot’ selection. In fact, the bases will be
constructed on equally spaced knots which have an additional computational
advantage. The approach is through utilizing the obtained irregular ‘knot’
distribution to change the topology of spline spaces. The two methods of doing
this will be considered: the first one is based on changing the topology of
the state space (the space of values of splines) and the second one is based
on transforming the domain space of splines.

The note is organized as follows. In the next section, we briefly recap basic
facts about the orthonormal basis of spline called splinets, which were
recently introduced in [5]. There, we also define two dimensional orthonormal
tensor splinets. This is followed by a section on the proposed method of the
knot selection in 2D. The two methods considered in this note are described in
the two subsequent sections. First, the state-space transformation method is
discussed in the 1D case and then explain the analogous approach to the 2D
case. A similar scheme is followed in the section on the domain transformation
method. A description of empirical studies that are designed to test the
methods is concluded in the note.

##  2 Tensor splinets

The splinets are efficient orthonormal bases of splines discussed in full
detail in [5] and implemented in the $R$-package Splinets, [6]. For the
purpose of this note, it is enough to know that a splinet is obtained by an
efficient dyadic algorithm from $B$-splines spanned on a given set of knots.
Two examples are shown in Figure 2, where the effect of different knot
distributions is clearly seen.

Figure 2:  Examples of splinets of the third order in the graphical
representation of the dyadic pyramid. In the figure, we see that placement of
knots changes the ability to model local variability in the data.

In the package, the focus is on one-dimensional domains. However, the tensor
product bases are simple enough that the package can also be utilized for the
analysis of functional data defined in a two-dimensional space.

Figure 3:  Top: Two $B$-splines (left) and corresponding $OB$-splines (right).
Bottom: The corresponding tensor spline $B_{\boldsymbol{\xi}i}^{k}\otimes
B_{\boldsymbol{\xi},i}^{l}$ (left) and the orthogonal tensor spline
$OB_{\boldsymbol{\xi}i}^{k}\otimes OB_{\boldsymbol{\xi},i}^{l}$(right).

Let us consider two knot selections $a=\xi_{0}<\dots<\xi_{n+1}=b$ and
$c=\zeta_{0}<\dots<\zeta_{m+1}=d$ on interval $[a,b]$ and $[c,d]$,
respectively. Further, let us consider two sets of $B$-splines
$\\{B_{\boldsymbol{\xi}i}^{k}(x)\\}$, $i=1,\dots,n-k+1$, $x\in[a,b]$ and
$\\{B_{\boldsymbol{\zeta}j}^{l}(y)\\}$, $j=1,\dots,m-l+1$, $y\in[c,d]$ of
orders $k$ and $l$, respectively, and spanned over the respective set of
knots. The linear functional spaces spanned on these sets of $B$-splines are
denoted by $\mathcal{H}^{k}_{\boldsymbol{\xi}}$ and
$\mathcal{H}^{l}_{\boldsymbol{\zeta}}$, respectively. Then the tensor product
$B$-spline basis is defined as

$\displaystyle\left(B_{\boldsymbol{\xi}i}^{k}\otimes
B_{\boldsymbol{\zeta}j}^{l}\right)(x,y)=B_{\boldsymbol{\xi}i}^{k}(x)B_{\boldsymbol{\zeta}j}^{l}(y),\,(x,y)\in[a,b]\times[c,d].$
(1)

Let
$\mathcal{H}^{k}_{\boldsymbol{\xi}}\otimes\mathcal{H}^{l}_{\boldsymbol{\zeta}}$
be the $(n-k+1)(m-l+1)$-dimensional linear functional space on
$[a,b]\times[c,d]$ spanned by the tensor product basis.

######  Theorem 2.1.

Let $f$ be square-integrable function on $[a,b]\times[c,d]$ and the splinets
$\\{OB_{\boldsymbol{\xi}i}^{k}\\}_{i=1}^{n-k+1}$ and
$\\{OB_{\boldsymbol{\xi},i}^{l}\\}_{j=1}^{m-l+1}$ are the orthogonormal
splines bases obtained by dyadic orthogonalization of
$B_{\boldsymbol{\xi}i}^{k}$ and $B_{\boldsymbol{\xi}j}^{l}$, respectively.
Then

$\left\\{OB_{\boldsymbol{\xi}i}^{k}\otimes
OB_{\boldsymbol{\zeta}j}^{l}\right\\}_{i,j=1}^{n-k+1,m-l+1}$

also constitute an orthonormal basis in
$\mathcal{H}^{k}_{\boldsymbol{\xi}}\otimes\mathcal{H}^{l}_{\boldsymbol{\zeta}}$
and the projection $\hat{f}$ of $f$ to this space is given by

$\hat{f}=\sum_{i=1}^{n-k+1}\sum_{j=1}^{m-l+1}\alpha_{ij}\cdot
OB_{\boldsymbol{\xi}i}^{k}\otimes OB_{\boldsymbol{\zeta}j}^{l},$

where
$\alpha_{ij}=\int_{a}^{b}OB_{\boldsymbol{\xi}i}^{k}(x)\int_{c}^{d}OB_{\boldsymbol{\zeta}j}^{l}(y)f(x,y)~{}dy~{}dx$.

Moreover, for each $f$ satisfying zero boundary conditions of the order $l$
and $m$ along $[a,b]\times\\{0\\}$ and $\\{0\\}\times[c,d]$, respectively, any
$k$ and $l$, each $\epsilon>0$ there exist knot selections
$a=\xi_{0}<\dots<\xi_{n+1}=b$ and $c=\zeta_{0}<\dots<\zeta_{m+1}=d$ such that

$\left(\int_{a}^{b}\int_{c}^{d}\left|f(x,y)-\hat{f}(x,y)\right|^{2}~{}dx~{}dy\right)^{1/2}\leq\epsilon.$

In Figure 3, we present abscissa and ordinate 1D-basis elements at the top and
$B_{\boldsymbol{\xi}i}^{k}\otimes B_{\boldsymbol{\zeta}j}^{l}$ (left) and
$OB_{\boldsymbol{\xi}i}^{k}\otimes OB_{\boldsymbol{\zeta}j}^{l}$ (right) at
the bottom.

##  3 ‘Knot’ selection in 2D

The algorithm that has been used for the 1D knot selection is applicable to
functional data that are defined on the 2D plane. Indeed, the method is
equivalent to building a regression binary tree where the so-called split
points (the nods at the tree) are found based on the best mean square error
(MSE) approximation of a regression function by functions that are constant
over a 2D grid, see [7]. The centers of rectangles of the 2D grid define
‘knots’ as shown in Figure 4 (left). The ‘knots’ for the MNIST fashion images
and their gradient images are presented in Figure 4 (middle, right), see [8]
for more details on this data set.

Figure 4:  (Left): Illustration of how the ‘knots’ in 2D are obtained by
taking the midpoints of 2D rectangular grid cells obtained by the binary
regression tree MSE fit to the functional data. (Middle) The original image
and the selected ‘knots’, (Right) The gradient image with ‘knots’.

We are using quotation marks when referring to obtained points since they are
not convenient to build direct spline approximation in a similar manner as in
the one-dimensional case. The reason for this is that there are no convenient
spline bases spanned on irregularly distributed ‘knots’ in such domains. One
could consider triangulation but defining splines on a triangular grid is
complicated and the resulting spline bases are difficult to orthogonalize. The
easiest to use and most popular 2D space of splines are based on the tensor
product bases, see, for example, [9]. However, the tensor bases are defined
over lattice grids, and thus the 2D grids shown in Figure 4 (Bottom) cannot be
directly used for constructing such bases.

##  4 The state space transformation

The problem of irregularly spread ‘knots’ occurs only for domains of a
dimension higher than one. Our proposed methods account for ‘knots’
distribution not directly in the construction of spline bases and can be
applied in 1D as well. In this case, the irregular knots can be used directly
in building splines thus one can benchmark the proposed methods against this
direct irregular knot selection.

The main component of both methods is a continuous approximation of the
irregular ‘knot’ distribution. We consider ‘knots’ $\boldsymbol{\xi}_{j}$,
$j=1,\dots,N$ that are irregularly distributed. One can simply apply kernel-
driven estimation-like methods known in non-parametric statistics by treating
the ‘knots’ as iid data drawn from a certain distribution. Since ‘knots’ are
not truly iid data, those methods should be viewed as ad hoc smoothing
methods. More suitably, one can use Splinet-package to project discretized
data to the space of splines as presented in Theorem 2.1. From now on, $g$
represents the density of the ‘knot’ distribution.

Small values of $g$ indicate regions where the data are less informational.
This motivates us to introduce a $g$-weighted topology in $L^{2}_{g}=\\{h:\int
h^{2}g<\infty\\}$

$\langle h,h\rangle_{g}=\int
h^{2}(x)g(x)dx=\int\left(h(x)\sqrt{g(x)}\right)^{2}dx=\langle
h\sqrt{g},h\sqrt{g}\rangle,\ \ \ \forall h\in L^{2}_{g}.$

We can transfer the functional data $X_{1},X_{2},\dots,X_{n}$ to

$\tilde{X_{i}}(x)=X_{i}(x)\sqrt{g(x)}.$ (2)

One can use $\tilde{X_{1}}(x),\tilde{X_{2}}(x),\dots,\tilde{X_{n}}(x)$ instead
of the original data and perform any functional analysis with one important
modification. Namely, since the knot distribution is already accounted for,
one does not need to use splines built on irregularly distributed knots. This
is crucial in the 2D case since the equally spaced knots constitute a lattice
on which tensor splines are well defined.

######  Remark 4.1.

The sparsity of ‘knots’ in some regions may lead to zero (or nearly zero)
values of $g$. This may lead to a subdomain of the original
$[a,b]\times[c,d]$. This subdomain could constitute the support of the
functional data. To benefit from that one would like to have the possibility
to account for a reduced support set. This feature is actually implemented in
the current version of Splinets-package but will be extended to the two-
dimensional case in order to utilize sparsity of the images.

##  5 The domain space transformation

Here, we also consider the ‘knots’ distribution but instead of the density
$g$, we utilize the distribution function
$G(x,y)=\int_{a}^{x}\int_{c}^{y}g(u,v)~{}du~{}dv$. For that $g$ approximated
in the previous method can be used. Alternatively, one can direct approach to
the problem by looking at the ‘empirical’ equivalent

$\hat{G}(x,y)=\frac{\\#\\{{\boldsymbol{\xi}}\cap[a,x]\times[c,y]\\}}{N}.$

Then $\hat{G}$ can be used to construct domain mapping $\mathbf{G}$ from
$[a,b]\times[c,d]$ to $[0,1]\times[0,1]$ such that
$\mathbf{G}(\boldsymbol{\xi})$ are uniformly distributed.

To simplify the discussion in this note we consider only the one-dimensional
case. Let $\nu_{i}={G}_{X}(\xi_{i})$, yielding approximately uniformly
distributed $\nu_{1},\ldots,\nu_{n}$. Thus our functional data point, say $h$,
needs to be transformed to $\tilde{h}=h\circ G^{-1}$. For these domain-
transformed data, we have

$\displaystyle\langle\tilde{h},\tilde{k}\rangle$
$\displaystyle=\int_{0}^{1}h(G^{-1}(u))k(G^{-1}(u))~{}du=\int_{0}^{1}h(x)k(x){g(x)}~{}du=\langle
h,k\rangle_{g}$

We conclude that while the two methods transform data in different ways they
are topologically equivalent.

######  Theorem 5.1.

Let $g$ be the density distribution on $[0,1]$, $G$ be the corresponding cdf.
Moreover, let $h_{g}=\sqrt{g}h$ for $h\in L_{g}^{2}$ and $\tilde{h}_{G}=h\circ
G^{-1}$, then

$\langle h_{g},k_{g}\rangle=\langle\tilde{h}_{G},\tilde{k}_{G}\rangle.$

Figure 5: Left: Approximated density of knots for the data-driven approach.
The knots were selected using the DDK package. Right: One-dimensional Hilbert
curve representations of one sample from the T-shirt class from the MNIST
data.

##  6 Data analysis

To assess the efficiency we utilize the fashion MNIST image dataset. For
simplicity, we focus on one class of this dataset, the T-shirt class. The
performance was compared with the conventional equidistant and the full data-
driven knot selection methodologies detailed in [1]. Since we focus in this
note on the 1D analysis (leaving a fully 2D discussion for a full paper), we
looked at the images only as one-dimensional vectors. However, to better
account for spatial dependencies in the images, we used the Hilbert curve to
order the pixels in this vector, for more details see [10].

Our experiment was organized into three distinct workflows. In the first
workflow (W1), we employ the traditional strategy of equidistant knot
selection. Subsequently, in the second workflow (W2), we leveraged the Data-
Driven Knot (DDK) selection algorithm detailed in the aforementioned citation
to select the knots. Finally, in the third workflow (W3), which encapsulates
our proposed methodology, we utilized the knots selected from the DDK
algorithm to generate the corresponding density function using the KernSmooth
R package. The image data was then weighted by this density function,
following the process outlined in (2).

In all three workflows, following the selection of the knots, we employ the
Splinets R package to project the data onto third-order splines with the
selected knots. Subsequently, we perform a spectral decomposition of the data
by estimating the eigenvalues, denoted $\lambda_{i}$, and the corresponding
eigenfunctions. Next, we compute the explained variance for each of the three
workflows, labeled $EV_{1}$, $EV_{2}$, and $EV_{3}$. To achieve this, we first
perform PCA on our original data set to obtain a set of eigenvalues
$\Lambda_{k}$, $k=1,\ldots,784$.

Figure 6: Explaind variance for the three workflows

For each workflow, we obtained a set of eigenvalues, which are derived from
the projection of the data to the corresponding workflow subspace. That is for
$i=1,2,3$ we get $\lambda_{j}^{(i)}$, $j=1,\ldots,n$, where $n$ is a pre-
specified number of eigenvalues that we want the explained variance with
respect to. The explained variance $EV_{i}$ for each workflow $i$ is then
calculated using

$EV_{i}=\frac{\sum_{j=1}^{n}\lambda_{j}^{(i)}}{\sum_{k=1}^{784}\Lambda_{k}.}$

This equation quantifies the proportion of the total variance in the original
data that is captured by $n$ eigenvalues in the subspace resulting from each
workflow. In Figure 6, we observe the $W2$ and $W3$ provide better-explained
variance compared to $W1$. For example, for $n=10$ we have
$EV_{1}=0.547,EV_{2}=0.635,EV_{3}=0.642$.

## References

  * [1] Basna, R., Nassar, H., Podgórski, K.: Data driven orthogonal basis selection for functional data analysis.  Journal of Multivariate Analysis 189, 104,868 (2022)
  * [2] Basna, R., Nassar, H., Podgórski, K.: R-Package: DKK, Version: 0.1.0, Orthonormal Basis Selection using Machine Learning (2021).  URL https://github.com/ranibasna/ddk
  * [3] Meurens, M.: Wine – data set.  Online.  URL https://rdrr.io/cran/cggd/man/Wine.html
  * [4] Podgórski, K.: Splinets – splines through the taylor expansion, their support sets and orthogonal bases (2021)
  * [5] Liu, X., Nassar, H., Podgórski, K.: Dyadic diagonalization of positive definite band matrices and efficient b-spline orthogonalization.  Journal of Computational and Applied Mathematics 414, 114,444 (2022).  https://doi.org/10.1016/j.cam.2022.114444.  URL https://www.sciencedirect.com/science/article/pii/S0377042722002102
  * [6] Liu, X., Nassar, H., Podgórski, K.: Splinets: Functional Data Analysis using Splines and Orthogonal Spline Bases.  CRAN.R-project.org (2023).  R package version 1.5.0
  * [7] Hastie, T., Tibshirani, R., Friedman, J.: The elements of statistical learning: data mining, inference and prediction, 2 edn.  Springer (2009).  URL http://www-stat.stanford.edu/~tibs/ElemStatLearn/
  * [8] Fashion MNIST - data set.  Online, Zalando Research.  URL https://www.kaggle.com/datasets/zalando-research/fashionmnist?resource=download
  * [9] Eilers, P., Marx, B.: Practical smoothing. The joys of $P$-splines.  Cambridge University Press, Cambridge (2021)
  * [10] Bader, M.: Space-filling curves: an introduction with applications in scientific computing, vol. 9.  Springer Science & Business Media (2012)