aoraoraorist attrattrattributive

# Enhancements to the BOUN Treebank Reflecting the Agglutinative Nature of
Turkish

Büşra Marşan1
<EMAIL_ADDRESS>Salih Furkan Akkurt2
<EMAIL_ADDRESS>Muhammet Şen2
<EMAIL_ADDRESS>Merve Gürbüz2
<EMAIL_ADDRESS>Onur Güngör2
<EMAIL_ADDRESS>Şaziye Betül Özateş2
<EMAIL_ADDRESS>Suzan Üsküdarlı2
<EMAIL_ADDRESS>Arzucan Özgür2
<EMAIL_ADDRESS>Tunga Güngör2
<EMAIL_ADDRESS>Balkız Öztürk1
<EMAIL_ADDRESS>
1\. Linguistics Department, Boğaziçi University
2\. Computer Engineering Department, Boğaziçi University

###### Abstract

In this study, we aim to offer linguistically motivated solutions to resolve
the issues of the lack of representation of null morphemes, highly productive
derivational processes, and syncretic morphemes of Turkish in the BOUN
Treebank without diverging from the Universal Dependencies framework. In order
to tackle these issues, new annotation conventions were introduced by
splitting certain lemmas and employing the MISC (miscellaneous) tab in the UD
framework to denote derivation. Representational capabilities of the re-
annotated treebank were tested on a LSTM-based dependency parser and an
updated version of the BoAT Tool is introduced.

_Keywords_ Universal Dependencies $\cdot$ Turkish $\cdot$ morphological
analysis $\cdot$ dependency annotation $\cdot$ dependency parsing

##  1 Introduction

Following the dependency grammar framework first proposed by Tesniére [20],
dependency trees illustrate how sentence elements relate to one another
through head and dependent relations. Universal
Dependencies111https://universaldependencies.org (UD) is an international
cooperative treebank project based on the dependency grammar framework and it
aims to offer a standardized and comprehensive dependency treebank collection
covering 121 languages.

With the addition of new UD treebanks, Turkish does not qualify as a low
resource language anymore. With a total of 733,000 tokens, it is the 12th
largest UD treebank in the UD repository. Although the coverage of the
treebanks plays an essential role in improving the performance of natural
language processing (NLP) systems [8], their ability to correctly and
consistently illustrate the morphosyntactic features of the target language
should not be overlooked. As Vincze et al. [22]’s study shows, the better a
treebank’s ability to represent the morphology and syntax of the target
language, the better the performance of the NLP systems using that treebank as
a resource.

In this paper, we aim to abide by the linguistic framework set by Bedir et al.
[3] and offer an updated and comprehensive UD treebank for Turkish, the BOUN
Treebank, along with an improved UD annotation interface, the BoAT Tool, first
introduced in Türk et al. [21].

The decisions made in the re-annotation process of the BOUN Treebank aim to
offer solutions to the issues posed by the morphologically rich and complex
nature of Turkish: null morphemes are frequently employed, agglutinative
processes are heavily used to create new forms, and numerous morphemes like
copula and -ki are very syncretic. The main goal of this study is to
illustrate these phenomena without compromising the compliance with the UD
framework.

This paper is organized as follows. Previous attempts at creating dependency
treebanks in Turkish are laid out in Section 2. Annotation changes made in the
current version of the BOUN Treebank and their linguistic justification are
discussed in Section 3. Statistics regarding the changes made to the previous
version are stated in Section 4. Improvements made to the BoAT Annotation Tool
are explained in Section 5. Finally, the performance of the parser trained
using the current version of the treebank is reviewed in Section 6.

##  2 Dependency Treebanks in Turkish

Shortly after the first dependency treebank for Turkish was presented by
Atalay et al. [2], Eryiğit and Pamay [6] offered a smaller dependency treebank
consisting of 300 sentences as part of the CoNLL 2007 Shared Task: MST
Treebank. 13 years after its publication, Sulubacak et al. [18] re-annotated
this dataset, converted it to the UD framework, and published the updated
dataset as IMST-UD. Çöltekin’s The Grammar Book Treebank (GB) [4] which
consists of 2,803 sentences extracted by a reference book on Turkish grammar
by Göksel and Kerslake [9] marks the very first effort in creating the first
UD-style Turkish treebank. Another pioneer in Turkish dependency treebanks is
IWT-UD as it is the first Turkish dependency treebank that covers informal
texts. IWT-UD was introduced by Sulubacak and Eryiğit [17] three years after a
constituency-style treebank, IWT, was presented by Pamay et al. [15].

Tourism is one of the two domain-specific dependency treebanks in Turkish and
consists of hotel and restaurant reviews. The other one is ATIS treebank that
covers the Turkish translation of English ATIS (Airline Travel Information
System) corpus [11].

Other Turkish UD-style treebanks include Kenet UD Treebank, Penn Treebank, and
FrameNet treebank. Penn Treebank consists of Turkish translations of English
Penn Treebank [19] while FrameNet consists of 2,700 sentences from the Turkish
FrameNet database [13].

With 9,761 sentences and 121,214 tokens randomly selected from Turkish
National Corpus (TNC) [1], the BOUN Treebank is one of the largest UD-style
treebanks in Turkish. Covering five different registers (broadsheet national
newspapers, biographical texts, essays, popular culture articles, and
instructional texts), it offers word order and sentence length variance in
addition to linguistically motivated dependency annotations (see Section 3).

##  3 Improving BOUN Treebank

The first step of the previous annotation process of the BOUN Treebank (see
[21] for a detailed discussion) was parsing the raw text to create CoNLL-U
files using Kanerva et al.’s [12] pipeline tool. During this parsing process,
UPOS tags and certain morphological information were automatically annotated.
Then the dependency relations were manually annotated by two native speakers
of Turkish who are linguists. To ensure inter-annotator agreement, randomly
selected 1000 sentences were double annotated. Using Cohen’s Kappa measure,
inter-annotator agreement was calculated. Dependency label match score was
0.82, unlabelled attachment score was 0.83, and labelled attachment score was
0.75.

###  3.1 The Re-annotation Process: Overcoming the Challenges

For the re-annotation process, a team of linguists detected shortcomings and
problematic annotations of the previous version of the BOUN Treebank [21]. Two
major representation challenges were detected: derivation and null copula. The
following subsections discuss the strategies employed to overcome these
challenges.

Annotations were done by two linguists who are native speakers of Turkish. In
order to ensure inter-annotator agreement, 100 sentences were double
annotated. Unlabelled attachment and labelled attachment scores were
calculated using Cohen’s Kappa measure. Their respective values are 98.61 and
97.81.

####  3.1.1 Derivation

Having its focus on syntax, the UD framework falls short of representing
derivational processes. The official guide of UD argues that the final
derivational suffix in a word is opaque in the sense that it does not permit
access to the morphemes that come before it. Hence in a construction as the
one shown in Figure 1, numerous derivational and inflectional morphemes before
the last derivational morpheme (-ki) are lost.

Figure 1: Derivational and inflectional analysis of “görüşmelerindeki” (“those
who were/are at his/her meetings”)

As one of the primary concerns of this study was finding ways to illustrate
derivation processes like those in Figure 1 without diverging from the UD
framework, two strategies were employed for different cases:

  * •

For morphemes like -lI (“with”), df= function is introduced in the MISC tab.

  * •

Lemmas containing -ki morpheme are splitted.

These derivational morphemes imply that the host lemma and its modifier form a
syntactic unit together. As a result, the correct bracketing of the -lI
adjective, its modifier, and the constituent they modify together should look
like (3.1.1) instead of (3.1.1). Representing the derivational morphemes like
-lI and -sIz allows keeping this crucial syntactic
information.222Abbreviations: attr = attributive

[ [ [kahverengi] tüy-lü] kedi]
[brown fur- cat
“a cat with brown fur” [ [kahverengi tüy]-lü kedi]
[brown fur- cat
“a cat with brown fur”

With the new df= function proposed in this study, “tüylü” is not decomposed as
two lemmas: “tüy” (“fur”) and “lü” (“with”). Instead, it is left intact and
df=tüy (“derived from=fur”) function is annotated in the MISC tab.

Another challenge regarding the derivational processes is posed by -ki. There
are two different -ki morphemes in Turkish [10]. One is used as noun and the
other derives adjectives from nouns.

Statistically, the vast majority of inflectional suffixes in Turkish occur
after the derivational ones. However, both types of the -ki morphemes diverge
from this distribution as it can be observed in Figure 1, where the
inflectional suffix -ler precedes the derivational -ki. Hence allowing -ki to
“block access” [3]333It is stated in the UD guidelines that “the lemma does
not remove derivational morphology, so the lemma of [en] “organizations” is
“organization” not “organize” (nor “organ”).” See
https://universaldependencies.org/u/overview/morphology.html for the full
picture. to the morphemes that were attached before it reduces the
capabilities of the annotation to represent a great deal of morphosyntactic
information. With the aim of offering a solution to this issue, a decision to
split lemmas that contain either type of -ki was made.444Refer to the Appendix
to compare previous and updated annotation schemes for both -ki morphemes.
Although splitting lemmas is not a common practice within the UD framework, we
believe that the theoretical motivation behind this decision justifies the
divergence from the framework. (For a detailed linguistic discussion of this
issue, please refer to Bedir et al. [3])

####  3.1.2 Null Morpheme

Languages such as Turkish, Russian, Arabic, and Coptic have null morphemes,
however, the UD framework does not officially support such phenomena. As a
result, independent strategies have emerged to represent null morphemes in
these languages. For example, Coptic avoids null subject nodes by using fused
forms[23], Marathi annotates the feature values of the null morphemes,
however, does not introduce any information (i.e. annotation) to indicate that
they are null morphemes[16]. Widely employing null morpheme for pluralization,
Arabic distinguishes between two types of annotation: Form-based and function-
based. In their UD-style dependency treebank annotation, Marton et al.[14]
follow a function-based annotation framework and annotate the feature values
of null morphemes.

In Turkish, the copula can surface in three different forms [9]: i-, -y-, and
$\emptyset$.After considering UD guidelines and particularities of Turkish
copula, it was decided to employ the MISC tab again by introducing two new
functions for the null copula: nullcop=3s (singular) and nullcop=3p (plural).
By following a function-based annotation schema, we were able to offer more
linguistically accurate annotations without diverging from the UD framework.

####  3.1.3 Copula

In Turkish, ol- copula has six distinct functions [3]: An intransitive verb
meaning “to be suitable/fit”, a transitive verb meaning “to become”, an
auxiliary verb in embedded sentences, an auxiliary verb following the
participle, a light verb forming complex verbal constructions (such as “sorun
olmak” (“to be/become an issue”)), and finally the existential predicate that
surfaces as “var” (“to exist”) and “yok” (“not to exist”). Yet the previous
annotation scheme of the BOUN Treebank made no distinction between these
different usages. To offer more accurate representations, certain annotation
changes were made regarding these functions.

###  3.2 Newly Introduced XPOS Tags and Dependency Relations

In an attempt to overcome the challenges thoroughly discussed in the previous
subsection, a set of new XPOS tags and dependency relations were introduced in
the updated version of the BOUN Treebank. A comprehensive list can be found in
Table 1.

Lemma | var | yok |  | ol- (after
---
participle)
| ol- (in
---
embedded
sentences)
| ol- (in light verb
---
constructions)
| ol- (as transitive
---
or intransitive
verb)
| -ki
---
(adjectivizer)
| -ki
---
(pronominal)
UPOS | NOUN | NOUN | AUX | AUX | VERB | VERB | PART | PRON
XPOS | Exist | Exist |  | Ptcp |  |  | Attr | Partic
Deprel | root | root | aux | cop | compound:lvc | root | dep:der |

Table 1: New dependency relations and XPOS tags proposed for the updated BOUN
Treebank.

##  4 Statistics

As part of the re-annotation process, 117,732 changes were made in the
following tabs: UPOS, XPOS, Deprel, MISC, and Features. The majority of
annotation changes targeted UPOS and XPOS tags (see Table 2). Since
morphological information was automatically annotated in the previous version
of the treebank by the parsing tool [12], refinements by the annotators were
required in order to ensure accuracy.

Field | UPOS | XPOS | Features | Deprel | MISC
---|---|---|---|---|---
Changes | 11,396 | 63,829 | 27,098 | 23,32 | 4,973
Table 2: Changes made in the re-annotation process.

The changes in UPOS values reflect the linguistics-based decisions made in the
re-annotation process. Previous version of BOUN Treebank made no distinction
between two -ki morphemes in Turkish. As a result, almost all -ki instances
were labeled as CConj (clausal conjunction). After deciding to make a
distinction between adjectivizer -ki and pronominal -ki, the UPOS tag of the
former was changed to Part.

Field | UPOS | XPOS
---|---|---
Change | Adj ->Noun | CConj ->Part | Noun ->Propn | Verb ->Ptcp | Verb ->Vnoun | ANum ->Indef
Count | 1,595 | 1,025 | 968 | 2,459 | 1,664 | 1,622

Table 3: Most frequent changes targeting the UPOS and XPOS tags

Due to the shortcomings of automatic morphological tagging, some proper nouns
were labeled as nouns. Re-annotation process targeted them as well: UPOS tags
of 986 proper nouns were changed. In addition, XPOS tags of verbal nouns and
participle forms were updated (see Table 3).

##  5 The BoAT Tool

The BoAT tool offered in Türk et al. [21] is a desktop application for
manually annotating sentences parsed by a dependency parser. In the scope of
the current work, in addition to the reannotation of the BOUN Treebank, we
enhanced the tool with additional functionalities. We will publish the tool as
open source on GitLab accompanied with a user manual.

###  5.1 Changes

BoAT is a desktop application, written in Python and based on Qt. The Qt
version has been incremented from 5 to 6. This resulted in a more modern-
looking user interface (UI). With ample feedback from the annotators who had
used the tool for the BOUN Treebank, some requested features and improvements
have been surfaced.

Clutter: Annotation table’s columns were being shown or hidden by checkboxes
above the table. These were removed and a textbox beside the other buttons has
been added to replace them. This textbox serves exactly the same purpose while
taking less space.

Autocompletion: Annotators use the tool for long hours at a time. Thus after a
while, mistakes tend to occur. Another requested feature, autocompletion of
the table fields, aims to prevent such mistakes. Many fields of the table have
predetermined sets of values they can take. By not allowing values outside
these sets and having a shorthand writing system whereby an annotator enters
only the start of a value and it gets filled automatically, we implemented
this much requested feature.

Saving as CoNLL-U: Another change was regarding saving of the CoNLL-U
documents. The initial tool saved every edit automatically, yet treebanks with
thousands of sentences tend to take time to save. This mishap seemed to slow
our annotators. We added a save button instead of the autosave feature.
Currently, a user uses the save button to edit the actual CoNLL-U file.

Shortcuts: The initial tool already had shortcut support. Going with the
focus-oriented approach, all the new tasks have keyboard shortcuts associated
with them as well. This alleviates the need to use a mouse while annotating
via keyboard.

Dependency graphs: The vertical dependency graph in the initial version was
replaced by the horizontal dependency graph of spaCy [7], displaCy, due to the
spatial concerns.

Validation: The validation script for annotations has been upgraded to the
latest version written by the UD framework, which has much more detailed
explanations for why a specific annotation is invalid.

##  6 Parser Performance

Within the scope of this study, an NLP task was conducted. BiLSTM-based
biaffine dependency parser proposed by Dozat and Manning [5] was trained using
the updated BOUN Treebank. The train set contained 7.803 sentences,
development set contained 982 sentences and test set contained 979 sentences.
Considering the average arc length and average token count (see Table 4 for
details) of each set, BOUN Treebank offers a well-balanced data set.

| Train | Development | Test | Entire Data
---|---|---|---|---
| Average
---
Arc Length
2.91 | 2.88 | 2.82 | 2.90
| Average
---
Token Count
12.83 | 12.42 | 12.36 | 12.74
| Number
---
of Sentences
7,803 | 982 | 979 | 9,761
Table 4: Sizes and specifications of train, development and test data sets

The previous version of the BOUN Treebank [21] yielded 77.36 unlabeled
attachment score (UAS) and 70.37 labeled attachment score (LAS). After being
trained on the new dataset, UAS is increased by 0.59 points to reach 77.96
while LAS showed a 0.10 decrease by dropping to 70.26 points. After the re-
annotation process, several faulty or disputed dependency relations were fixed
in the data, hence a rise in UAS was observed.

In order to better account for the particularities of Turkish morphosyntax,
four new dependency labels were introduced: dep:der for adjectivizer -ki,
obl:tmod for obliques that offer temporal information regarding the predicate,
advmod:emph for dA clitics, and compound:lvc for light verb constructions with
ol- copula. Moreover, new use cases for cop dependency label were offered to
represent different functions of the ol- copula. These changes in the
annotation framework added four new dependency types, thus the total class
number was increased.

A sum of 1,032 dep:der, 894 obl:tmod, 1,860 advmod:emph, and 1,545
compound:lvc tags were added. Newly introduced tags and classes introduced
added complexity, hence they might be the reason behind the slight decrease in
the LAS results.

The main annotation changes made as part of this study were focused on
morphology yet BiLSTM-based biaffine dependency parser proposed by Dozat and
Manning [5] doesn’t refer to the morphological information. In fact, it almost
completely ignores the morphology features while parsing. Hence the
significance of the improvements offered by this study can be better gauged by
a morphology-aware parser or a morphological analysis based downstream task.

##  7 Conclusion and Further Research

The UD framework highly emphasizes syntactic relations and aims to offer a
universal foundation to represent typologically different languages in a
uniform way. While doing so, certain particularities of these languages tend
to get lost in the annotation process in an attempt to abide by the UD
convention. The aim of this study is to offer linguistically sound solutions
to illustrate syntactically relevant morphological features of Turkish such as
null morpheme realizations and derivational processes without diverging
significantly from the UD framework.

In order to test the morphological capabilities of the re-annotated BOUN
Treebank, a parser that refers to morphological information can be implemented
in further research.

acknowledgments This work was supported by Boğaziçi University Research Fund
Grant Number 16909. TUBAGEBIP Award of the Turkish Science Academy (to A.O.)
is gratefully acknowledged.

## References

  * Aksan et al. [2012] Y. Aksan, M. Aksan, A. Koltuksuz, T. Sezer, Ü. Mersinli, U. U. Demirhan, H. Yılmazer, G. Atasoy, S. Öz, İ. Yıldız, et al.  Construction of the turkish national corpus (tnc).  In _Proceedings of the Eighth International Conference on Language Resources and Evaluation (LREC’12)_ , pages 3223–3227, 2012.
  * Atalay et al. [2003] N. B. Atalay, K. Oflazer, and B. Say.  The annotation process in the turkish treebank.  In _Proceedings of 4th International Workshop on Linguistically Interpreted Corpora (LINC-03) at EACL 2003_ , 2003.
  * Bedir et al. [2021] T. Bedir, K. Şahin, O. Güngör, S. Uskudarli, A. Özgür, T. Güngör, and B. Ö. Başaran.  Overcoming the challenges in morphological annotation of turkish in universal dependencies framework.  In _Proceedings of The Joint 15th Linguistic Annotation Workshop (LAW) and 3rd Designing Meaning Representations (DMR) Workshop_ , pages 112–122, 2021.
  * Cöltekin [2015] C. Cöltekin.  A grammar-book treebank of turkish.  In _Proceedings of the 14th workshop on Treebanks and Linguistic Theories (TLT 14)_ , pages 35–49, 2015.
  * Dozat et al. [2017] T. Dozat, P. Qi, and C. D. Manning.  Stanford’s graph-based neural dependency parser at the conll 2017 shared task.  In _Proceedings of the CoNLL 2017 shared task: Multilingual parsing from raw text to universal dependencies_ , pages 20–30, 2017.
  * Eryiğit and PAMAY [2007] G. Eryiğit and T. PAMAY.  Itu validation set for metu-sabancı turkish treebank.  _Türkiye Bilişim Vakfı Bilgisayar Bilimleri ve Mühendisliği Dergisi_ , 7(1):31–37, 2007.
  * Explosion [2022] Explosion.  spaCy, 2022.  URL https://spacy.io/.  [Online; last accessed 27 March 2022].
  * Foth et al. [2014] K. Foth, A. Köhn, N. Beuck, and W. Menzel.  Because size does matter: The hamburg dependency treebank.  2014\.
  * Göksel and Kerslake [2004] A. Göksel and C. Kerslake.  _Turkish: A comprehensive grammar_.  Routledge, 2004.
  * Hankamer [2004] J. Hankamer.  Why there are two–ki’s in turkish.  _Current research in Turkish linguistics, Eastern Mediterranean University_ , 2004.
  * Hemphill et al. [1990] C. T. Hemphill, J. J. Godfrey, and G. R. Doddington.  The atis spoken language systems pilot corpus.  In _Speech and Natural Language: Proceedings of a Workshop Held at Hidden Valley, Pennsylvania, June 24-27, 1990_ , 1990.
  * Kanerva et al. [2018] J. Kanerva, F. Ginter, N. Miekka, A. Leino, and T. Salakoski.  Turku neural parser pipeline: An end-to-end system for the conll 2018 shared task.  In _Proceedings of the CoNLL 2018 Shared Task: Multilingual parsing from raw text to universal dependencies_ , pages 133–142, 2018.
  * Marşan et al. [2021] B. Marşan, N. Kara, M. Özçelik, B. N. Arıcan, N. Cesur, A. Kuzgun, E. Sanıyar, O. Kuyrukçu, and O. T. Yıldız.  Building the turkish framenet.  In _Proceedings of the 11th Global Wordnet Conference_ , pages 118–125, 2021.
  * Marton et al. [2011] Y. Marton, N. Habash, and O. Rambow.  Improving arabic dependency parsing with form-based and functional morphological features.  In _Proceedings of the 49th annual meeting of the association for computational linguistics: Human language technologies_ , pages 1586–1596, 2011.
  * Pamay et al. [2015] T. Pamay, U. Sulubacak, D. Torunoğlu-Selamet, and G. Eryiğit.  The annotation process of the itu web treebank.  In _Proceedings of the 9th Linguistic Annotation Workshop_ , pages 95–101, 2015.
  * Ravishankar [2017] V. Ravishankar.  A universal dependencies treebank for marathi.  In _Proceedings of the 16th International Workshop on Treebanks and Linguistic Theories_ , pages 190–200, 2017.
  * Sulubacak and Eryiğit [2018] U. Sulubacak and G. Eryiğit.  Implementing universal dependency, morphology, and multiword expression annotation standards for turkish language processing.  _Turkish Journal of Electrical Engineering & Computer Sciences_, 26(3):1662–1672, 2018.
  * Sulubacak et al. [2016] U. Sulubacak, M. Gökırmak, F. Tyers, Ç. Çöltekin, J. Nivre, G. Eryiğit, et al.  Universal dependencies for turkish.  In _Proceedings of COLING 2016, the 26th International Conference on Computational Linguistics: Technical Papers_. The Association for Computational Linguistics, 2016.
  * Taylor et al. [2003] A. Taylor, M. Marcus, and B. Santorini.  The penn treebank: an overview.  _Treebanks_ , pages 5–22, 2003.
  * Tesnière [1959] L. Tesnière.  Eléments de syntaxe structurale, paris: Klinck sieck.  _German translation:(1980) Grundzüge der strukturalen Syntax Stuttgart: Klett-Cotta_ , 1959.
  * Türk et al. [2021] U. Türk, F. Atmaca, Ş. B. Özateş, G. Berk, S. T. Bedir, A. Köksal, B. Ö. Başaran, T. Güngör, and A. Özgür.  Resources for turkish dependency parsing: Introducing the boun treebank and the boat annotation tool.  _Language Resources and Evaluation_ , pages 1–49, 2021.
  * Vincze et al. [2017] V. Vincze, K. I. Simkó, Z. Szántó, and R. Farkas.  Universal dependencies and morphology for hungarian-and on the price of universality.  Association for Computational Linguistics, 2017.
  * Zeldes and Abrams [2018] A. Zeldes and M. Abrams.  The coptic universal dependency treebank.  In _Proceedings of the Second Workshop on Universal Dependencies (UDW 2018)_ , pages 192–201, 2018.

##  Appendix A A comparison of previous and updated annotation schemes

Token | Form | Lemma | UPOS | XPOS | Features | Head | Deprel | MISC
---|---|---|---|---|---|---|---|---
1 | başındaki | başındaki | ADJ | Adj |  | 2 | amod |
3 | şapkası | şapka | NOUN |  |  Case=Nom—Number=Sing—Number[psor]= Sing—Person=3—Person[psor]=3 | 0 | root |

Table 5: Previous annotation scheme for “başındaki şapkası” (“the hat on
his/her head”)

Token | Form | Lemma | UPOS | XPOS | Features | Head | Deprel | MISC
---|---|---|---|---|---|---|---|---
1-2 | başındaki | başındaki |  |  |  |  |  |
1 | başında | baş | NOUN |  |  Case=Loc—Number=Sing—Number[psor]= Sing—Person=3—Person[psor]=3 | 3 | nmod |
2 | ki | ki | PART | Attr |  | 1 | dep:der |
3 | şapkası | şapka | NOUN |  |  Case=Nom—Number=Sing—Number[psor]= Sing—Person=3—Person[psor]=3 | 0 | root |

Table 6: Updated annotation scheme for “başındaki şapkası” (“the hat on
his/her head”)

Token | Form | Lemma | UPOS | XPOS | Features | Head | Deprel | MISC
---|---|---|---|---|---|---|---|---
1 | seninki | seninki | NOUN |  | Case=Nom—Number=Sing | 0 | root |

Table 7: Previous annotation scheme for “seninki” (“that of yours”)

Token | Form | Lemma | UPOS | XPOS | Features | Head | Deprel | MISC
---|---|---|---|---|---|---|---|---
1-2 | seninki | seninki |  |  |  |  |  |
1 | senin | sen | PRON | PERS | Case=Gen—Number=Sing—Person=2—PronType=Prs | 2 | nmod:poss |
2 | ki | ki | PRON | Partic | Case=Nom—Number=Sing | 0 | root |

Table 8: Updated annotation scheme for “seninki” (“that of yours”)


# Scalable mixed-domain Gaussian processes

Juho Timonen  Aalto University, Department of Computer Science  Harri
Lähdesmäki  Aalto University, Department of Computer Science

###### Abstract

Gaussian process (GP) models that combine both categorical and continuous
input variables have found use e.g. in longitudinal data analysis and computer
experiments. However, standard inference for these models has the typical
cubic scaling, and common scalable approximation schemes for GPs cannot be
applied since the covariance function is non-continuous. In this work, we
derive a basis function approximation scheme for mixed-domain covariance
functions, which scales linearly with respect to the number of observations
and total number of basis functions. The proposed approach is naturally
applicable to Bayesian GP regression with arbitrary observation models. We
demonstrate the approach in a longitudinal data modelling context and show
that it approximates the exact GP model accurately, requiring only a fraction
of the runtime compared to fitting the corresponding exact model.

##  1 INTRODUCTION

Gaussian processes (GPs) offer a flexible nonparametric way of modeling
unknown functions. While Gaussian process regression and classification are
commonly used in problems where the domain of the unknown function is
continuous, recent work has seen use of GP models also in mixed domains, where
some of the input variables are categorical or discrete and some are
continuous. Applications of mixed-domain GPs are found e.g. in Bayesian
optimization (Garrido-Merchán and Hernández-Lobato, 2020), computer
experiments (Zhang and Notz, 2015; Deng et al., 2017; Roustant et al., 2020;
Wang et al., 2021) and longitudinal data analysis (Cheng et al., 2019; Timonen
et al., 2021). For example in biomedical applications, the modeled function
often depends on categorical covariates, such as treatment vs. no treatment,
and accounting for such time-varying effects is essential. Since all commonly
used kernel functions (i.e. covariance functions) are defined for either
purely continuous or purely categorical input variables, kernels for mixed-
domain GPs are typically obtained by combining continuous and categorical
kernels through multiplication. Additional modeling flexibility can be
obtained by summing the product kernels as has been done in the context of GP
modeling for longitudinal data (Cheng et al., 2019; Timonen et al., 2021).

It is well known that exact GP regression has a theoretical complexity of
$\mathcal{O}(N^{3})$ and requires $\mathcal{O}(N^{2})$ memory, where $N$ is
the number of observations. This poses a computational problem which in
practice renders applications of exact GP regression infeasible for large
data. Various scalable approximation approaches for GPs have been proposed
(see e.g. (Liu et al., 2020) for a review). However, many popular approaches,
such as the inducing point (Snelson and Ghahramani, 2006; Titsias, 2009) and
kernel interpolation (Wilson and Nickisch, 2015) methods, can only be applied
directly if the kernel (i.e. covariance) function is continuous and
differentiable. In addition, they typically require the Gaussian observation
model, which is not appropriate for modeling for example discrete, categorical
or ordinal response variables. See Section 3 for a review of previous methods.

In this work, we present a scalable approximation scheme for mixed-domain GPs
and additive mixed-domain GPs, where the covariance structure depends on both
continuous and categorical variables. We extend the Hilbert space reduced-rank
approximation (Solin and Särkkä, 2019) for said additive mixed-domain GPs and,
making it applicable to e.g. analysis of large longitudinal data. The approach

  * •

scales linearly with respect to data set size

  * •

allows a wide variety of different categorical kernels

  * •

allows product kernels that consist of any number of continuous and
categorical kernels, as well as sums of such products

  * •

allows an arbitrary observation model and full Bayesian inference for the
model hyperparameters

To our knowledge, there are no existing approaches that satisfy these
conditions.

##  2 GAUSSIAN PROCESSES

###  2.1 Definition

A Gaussian process (GP) is a collection of random variables, any finite number
of which has a multivariate normal distribution (Rasmussen and Williams,
2006). A function $f$ is a GP

$f\sim\mathcal{GP}\left(m(\bm{x}),k(\bm{x},\bm{x}^{\prime})\right)$ (1)

with mean function $m(\bm{x})$ and kernel (or covariance) function
$k(\bm{x},\bm{x}^{\prime})$, if for any finite number of inputs
$\\{\bm{x}_{i}\in\mathcal{X}\\}_{i=1}^{N}$, the vector of function values
$\bm{f}=\left[f(\bm{x}_{1}),\ldots,f(\bm{x}_{N})\right]^{\top}$ follows a
multivariate normal distribution
$\bm{f}\sim\mathcal{N}\left(\bm{m},\bm{K}\right)$ with mean vector
$\bm{m}=\left[m(\bm{x}_{1}),\ldots,m(\bm{x}_{N})\right]^{\top}$ and $N\times
N$ covariance matrix $\bm{K}$ with entries
$\\{\bm{K}\\}_{ik}=k(\bm{x}_{i},\bm{x}_{k})$. The mean function is commonly
the constant zero function $m(\bm{x})=0$, and we have this convention
throughout the paper. The kernel function encodes information about the
covariance of function values at different points, and therefore affects the
model properties crucially.

###  2.2 Bayesian GP regression

In GP regression, the conditional distribution of response variable $y$ given
covariates $\bm{x}$ is modeled as some parametric distribution $p(y\mid
f(\bm{x}),\bm{\theta}_{\text{obs}})$, where $\bm{\theta}_{\text{obs}}$
represents possible parameters of the observation model. The function $f$ has
a zero-mean GP prior with covariance function
$k(\bm{x},\bm{x}^{\prime}\mid\bm{\theta}_{\text{GP}})$ that has
hyperparameters $\bm{\theta}_{\text{GP}}$. We focus on Bayesian GP modeling,
where in addition to the GP prior for $f$, we have a paramemeter prior
distribution $p(\bm{\theta})$ for
$\bm{\theta}=\left\\{\bm{\theta}_{\text{GP}},\bm{\theta}_{\text{obs}}\right\\}$.
Given $N$ observations $\mathcal{D}=\\{y_{n},\bm{x}_{n}\\}_{n=1}^{N}$ our goal
is to infer the posterior

$p\left(\bm{\theta},\bm{f}\mid\mathcal{D}\right)\propto
p\left(\bm{\theta},\bm{f}\right)\cdot p(\bm{y}\mid\bm{f},\bm{\theta}),$ (2)

where $\bm{y}=\left[y_{1},\ldots,y_{N}\right]^{\top}$. The part

$p\left(\bm{\theta},\bm{f}\right)=p\left(\bm{f}\mid\bm{\theta}\right)\cdot
p(\bm{\theta})$ (3)

is the prior and

$p(\bm{y}\mid\bm{f},\bm{\theta})=\prod_{n=1}^{N}p(y_{n}\mid
f(\bm{x}_{n}),\bm{\theta}_{\text{obs}})$ (4)

is the likelihood. This task often has to be done by sampling from
$p\left(\bm{\theta},\bm{f}\mid\mathcal{D}\right)$ using MCMC methods, which
requires evaluating the right-hand side of Eq. 2 (and possibly its gradient)
thousands of times. As the likelihood and parameter prior usually are
independent over each parameter and data point, they scale linearly and are
not a bottleneck. Instead, computing the GP prior density

$p\left(\bm{f}\mid\bm{\theta}\right)=\mathcal{N}\left(\bm{f}\mid\bm{0},\bm{K}\right),$
(5)

where the $N\times N$ matrix $\bm{K}$ has entries
$\\{\bm{K}\\}_{in}=k(\bm{x}_{i},\bm{x}_{n}\mid\bm{\theta}_{\text{GP}})$, is a
costly operation as evaluating the (log) density of the $N$-dimensional
multivariate normal distribution has generally $\mathcal{O}(N^{3})$ complexity
(see Suppl. Section 4). Furthermore, the matrix $\bm{K}$ takes
$\mathcal{O}(N^{2})$ memory.

An often exploited fact is that if the observation model (and therefore
likelihood) is Gaussian, $\bm{f}$ can be analytically marginalized and only
the marginal posterior $p\left(\bm{\theta}\mid\mathcal{D}\right)$ needs to be
sampled. This reduces the MCMC dimension by $N$ and likely improves sampling
efficiency, but one however needs to evaluate $p(\bm{y}\mid\bm{\theta})$ which
is again an $N$-dimensional multivariate Gaussian. The $\mathcal{O}(N^{3})$
complexity and $\mathcal{O}(N^{2})$ memory requirements therefore remain. In
this paper, we generally assume an arbitary observation model, and defer the
details of the Gaussian observation model until Suppl. Section 5.

###  2.3 Additive GP regression

In additive GP regression, the modeled function $f$ consists of additive
components so that $f=f^{(1)}+\ldots+f^{(J)}$, and each component
$j=1,\ldots,J$ has a GP prior

$f^{(j)}\sim\mathcal{GP}\left(0,k_{j}(\bm{x},\bm{x}^{\prime})\right),$ (6)

independently from other components. This means that the total GP prior is
$f\sim\mathcal{GP}\left(0,k(\bm{x},\bm{x}^{\prime})\right)$ with

$k(\bm{x},\bm{x}^{\prime})=\sum_{j=1}^{J}k_{j}(\bm{x},\bm{x}’).$ (7)

Furthermore, $\bm{f}^{(j)}\sim\mathcal{N}\left(\bm{0},\bm{K}^{(j)}\right)$ for
each vector

$\bm{f}^{(j)}=\left[f^{(j)}(\bm{x}_{1}),\ldots,f^{(j)}(\bm{x}_{N})\right]^{\top},$
(8)

where $j=1,\ldots,J$. The matrix $\bm{K}^{(j)}$ is defined so that its
elements are $\\{\bm{K}^{(j)}\\}_{ik}=k_{j}(\bm{x}_{i},\bm{x}_{k})$. This
means that the prior for $\bm{f}=\bm{f}^{(1)}+\ldots+\bm{f}^{(J)}$ is
$\bm{f}\sim\mathcal{N}\left(\bm{0},\bm{K}\right)$, where
$\bm{K}=\sum_{j=1}^{J}\bm{K}^{(j)}$.

Bayesian inference with MCMC for additive GP models requires sampling all
$\bm{f}^{(j)}$, meaning that adding one component increases the number of
parameters by $N$ (plus the possible additional kernel hyperparameters).
Moreover, the multivariate normal prior (Eq. 5) needs to be evaluated for each
component, adding the computational burden. In the case of Gaussian
likelihood, adding more components does not add any multivariate normal
density evaluations as it still needs to be done only for
$p(\bm{y}\mid\bm{\theta})$. Also the marginal posteriors of each
$\bm{f}^{(j)}$ are analytically available (see Suppl. Section 5).

###  2.4 Mixed-domain kernels for longitudinal data

Longitudinal data is common in biomedical, psychological, social and other
studies and consists of multiple measurements of several subjects at multiple
time points. In addition to time (often expressed as subject age), other
continuous covariates can be measured. Moreover, in addition to subject id,
other categorical covariates, such as treatment, sex or country can be
available. In the statistical methods literature, such data is commonly
modeled using generalized linear mixed effect models (Verbeke and Molenberghs,
2000). In recent work (Quintana et al., 2016; Cheng et al., 2019; Timonen et
al., 2021), longitudinal data has been modeled using additive GPs, where,
similar to commonly used linear models, each component $f^{(j)}$ is a function
of at most one categorical and one continuous variable. Each variable is
assigned a one-dimensional base kernel and for components that contain both a
continuous and categorical kernel, the kernel $k_{j}$ is their product. As the
total kernel $k$ is composed of the simpler categorical and continuous kernels
through multiplication and addition, it has a mixed domain.

These models have the very beneficial property that the effects of individual
covariates are interpretable. The marginal posterior distributions of each
component be studied to infer the marginal effect of different covariates. As
an example, if $k_{j}$ is just the exponentiated quadratic (EQ) kernel

$k_{\text{EQ}}(x,x^{\prime})=\alpha^{2}\exp\left(-\frac{(x-x^{\prime})^{2}}{2\ell^{2}}\right)$
(9)

and $x$ is age, the component $f_{j}$ can be interpreted as the shared effect
of age. On the other hand, if $k_{j}$ is the product kernel
$k_{\text{EQ}\times\text{ZS}}{\left((x,z),(x^{\prime},z^{\prime})\right)}=k_{\text{EQ}}(x,x^{\prime})k_{\text{ZS}}(z,z^{\prime})$
where

$k_{\text{ZS}}(z,z^{\prime})=\begin{cases}1\hskip 21.90874pt\text{ if
}z=z^{\prime}\\\ -\frac{1}{C-1}\text{ if }z\neq z^{\prime}\end{cases}$ (10)

is the zero-sum (ZS) kernel (Kaufman and Sain, 2010) for a categorical
variable $z$ that has $C>1$ categories, $f^{(j)}$ can be interpreted as the
category-specific effect of the continuous covariate $x$. This also has the
property that the effect sums to zero over categories at all values for $x$
(see Timonen et al. (2021) for proof), which helps in separating the category
effect from the shared effect, if a model has both.

Sometimes it is required to mask effects that are present for only a subset of
the individuals, such as the case individuals when the data also has a control
group. Using the kernel language, effect of component $j$ can be masked by
multiplying $k_{j}$ by a binary kernel $k_{\text{BIN}}(z,z^{\prime})$ which
returns 0 if either $z$ or $z^{\prime}$ takes a value in any of the masked
categories, and 1 otherwise.

##  3 RELATED RESEARCH

##### GPs and categorical inputs

A suggested approach to handle GPs with categorical covariates is to use a
one-hot encoding which turns a variable with $C$ categories into $C$ binary
variables, of which only one is on at a time, and then apply a continuous
kernel for them. Garrido-Merchán and Hernández-Lobato (2020) highlight that
the resulting covariance structure is problematic because it does not take
into account that only one of the binary variables can be one at a time. This
poorly motivated approach might have originated merely from the fact that
common GP software have lacked support for categorical kernels. We find it
more sensible to define kernels directly on categorical covariates, as that
way we can always impose the desired covariance structure.

Category-specific effects of a continuous covariate $x$ can be achieved also
by assigning independent GPs for the different categories. This way we have
only continuous kernel functions, and can possibly use scalable approaches
that are designed for them. This limited approach however cannot define any
additional covariance structure between the categories, such as the zero-sum
constraint (Eq. 10). The ZS kernel is a special case of compound symmetry
(CS), and for example Roustant et al. (2020) concluded that a CS covariance
structure was more justified than using only indenpendent GPs in their nuclear
engineering application.

Chung et al. (2020) developed a deep mixed-effect GP model that facilitates
individual-specific effects and scales as $\mathcal{O}(PT^{3})$, where $P$ is
the number of individuals and $T$ is the number of time points. Zhang et al.
(2020) handled categorical inputs by mapping them to a continuous latent space
and then using a continuous kernel. While this approach can detect interesting
covariance structures, it does not remove need to perform statistical modeling
with a predefined covariance structure as in Section 2.4. Another related non-
parametric way to model group effects is to use hierarchical generalized
additive models (Pedersen et al., 2019), as smoothing splines can be seen as a
special case of GP regression (Kimeldorf and Wahba, 1970).

##### Scalable GP approximations

A number of approximation methods exist that reduce the complexity of GP
regression to $\mathcal{O}(N\cdot M^{2})$, where $M$ controls the accuracy of
the approximation. Popular approaches rely on global sparse approximations
(Quiñonero Candela and Rasmussen, 2005) of the $N\times N$ covariance matrix
between all pairs of data points, using $M$ inducing points. The locations of
these inducing points are generally optimized using gradient-based continuous
optimization simultaneously with model hyperparameters, which cannot be
applied when the domain is not continuous. In Fortuin et al. (2021), the
inducing-point approach was studied in purely discrete domains and Cao et al.
(2015) presented an optimization algorithm that alternates between discrete
optimization of inducing points, and continuous optimization of the
hyperparameters. Disadvantages of this method are that it cannot find inducing
points outside of the training data, does not perform full Bayesian inference
for the hyperparameters, and assumes a Gaussian observation model.

A Hilbert space basis function approach for reduced-rank GP approximation in
continuous domains, on which this work is based, was proposed by Solin and
Särkkä (2019). Its use in the practice in the Bayesian setting was studied
more in Riutort-Mayol et al. (2020).

##  4 MIXED-DOMAIN COVARIANCE FUNCTION APPROXIMATION

###  4.1 Basic idea

We continue with the notation established in Section 2, and note that $\bm{x}$
denotes a general input that can consist of both continuous and categorical
dimensions. We consider approximations $\bm{f}\approx\tilde{\bm{f}}$ that
decompose the GP kernel function as

$k(\bm{x},\bm{x}^{\prime})\approx\tilde{k}(\bm{x},\bm{x}^{\prime})=\sum_{m=1}^{M}\psi_{m}(\bm{x})\psi_{m}(\bm{x}^{\prime}),$
(11)

where functions $\psi_{m}$ have to be designed so that the approximation is
accurate but easy to compute. This is useful in GP regression, because we get
a low-rank approximate decomposition for the kernel matrix
$\bm{K}\approx\tilde{\bm{K}}=\bm{\Psi}\bm{\Psi}^{\top}$, where $\bm{\Psi}$ is
the $N\times M$ matrix with elements
$\left[\bm{\Psi}\right]_{n,m}=\psi_{m}(\bm{x}_{n})$. Using this approximation,
we can write the approximate GP prior
$\tilde{\bm{f}}\sim\mathcal{N}\left(\bm{0},\bm{\Psi}\bm{\Psi}^{\top}\right)$
using $M$ parameters $\xi_{m}$ with independent standard normal priors,
connected to $\bm{f}$ through the reparametrization
$\tilde{\bm{f}}=\bm{\Psi}\bm{\xi}$, where
$\bm{\xi}=\left[\xi_{1},\ldots,\xi_{M}\right]^{\top}$. Evaluating the prior
density $p(\bm{\xi})=\prod_{m=1}^{M}\mathcal{N}\left(\xi_{m}\mid 0,1\right)$
has now only $\mathcal{O}(M)$ cost. After obtaining posterior draws of
$\bm{\xi}$, we can obtain posterior draws of $\tilde{\bm{f}}$ with
$\mathcal{O}(NM)$ cost, which comes from computing the matrix $\bm{\Psi}$. The
likelihood (Eq. 4) can then be evaluated one data point at a time and the
total complexity of the approach is only $\mathcal{O}(NM+M)$. Furthermore, the
memory requirement is reduced from $\mathcal{O}(N^{2})$ to $\mathcal{O}(NM)$,
since we only need to store $\bm{\Psi}$ and never compute $\bm{K}$. This is
the approach used throughout this paper, and the focus is on how to design the
functions $\psi_{m}$ for different kernel functions so that the approximation
is accurate with $M\ll N$.

###  4.2 Continuous isotropic covariance functions

A continuous stationary covariance function
$k(\bm{x},\bm{x}^{\prime}):\mathbb{R}^{D}\times\mathbb{R}^{D}\rightarrow\mathbb{R}$
depends only on the difference $\bm{r}=\bm{x}-\bm{x}^{\prime}$ and can
therefore be written as $k(\bm{r})$. Such covariance functions can be
approximated by methods that utilize the spectral density

$\displaystyle S(\omega)$ $\displaystyle=\int
k(\bm{r})e^{-i\omega^{\top}\bm{r}}\text{d}\bm{r}.$ (12)

If the covariance function is isotropic, meaning that it depends only on the
Euclidean norm $\lVert\bm{r}\rVert$, also $S(\omega)$ is isotropic and can be
written as $S\left(\lVert\Omega\rVert\right)$, i.e. as a function of one
variable. As shown in (Solin and Särkkä, 2019), an isotropic covariance
function can be approximated as

$\tilde{k}(\bm{x},\bm{x}^{\prime})=\sum_{b=1}^{B}s_{b}\phi_{b}(\bm{x})\phi_{b}(\bm{x}^{\prime}),$
(13)

where $s_{b}=S(\sqrt{\lambda_{b}})$ and $\phi_{b}$ and $\lambda_{b}$ are the
$B$ first eigenfunctions and eigenvalues of the Dirichlet boundary value
problem

$\begin{cases}\frac{\partial^{2}}{\partial\bm{x}^{2}}\phi_{b}(\bm{x})&=\lambda_{b}\phi_{b}(\bm{x}),\hskip
28.45274pt\bm{x}\in\Omega\\\ \phi_{b}(\bm{x})&=0,\hskip
59.75095pt\bm{x}\in\partial\Omega\end{cases}$ (14)

for a compact set $\Omega\subset\mathbb{R}^{D}$. We see that this
approximation has the same form as Eq. 11 with $M=B$ and
$\psi_{m}(\bm{x})=\sqrt{s_{m}}\phi_{m}(\bm{x})$. The spectral density has a
closed form for many kernels, and the domain $\Omega$ can be selected so that
the eigenvalues and eigenfunctions have one too. Functions $\psi_{m}(\bm{x})$
are therefore easy to evaluate and the computation strategy described in
Section 4.1 can then be used. As an example, when $D=1$ and $\Omega=[L,L]$
with $L>0$, we have

$\begin{cases}\phi_{b}(\bm{x})&=\frac{1}{\sqrt{L}}\sin\left(\frac{\pi
b(\bm{x}+L)}{2L}\right)\\\ \lambda_{b}&=\left(\frac{\pi
b}{2L}\right)^{2}\end{cases}.$ (15)

and it was proven in (Solin and Särkkä, 2019) that in this case
$\lim_{L\rightarrow\infty}\left[\lim_{B\rightarrow\infty}\tilde{k}(\bm{x},\bm{x}^{\prime})\right]=k(\bm{x},\bm{x}^{\prime})$
uniformly for any stationary $k$ that has a regular enough spectral density.
For example for the EQ kernel (Eq. 9)
$k_{\text{EQ}}(r)=\alpha^{2}\exp\left(-\frac{r^{2}}{2\ell^{2}}\right)$, the
spectral density is
$S_{\text{EQ}}(\omega)=\alpha^{2}\ell\sqrt{2\pi}\exp\left(-\frac{\ell^{2}\omega^{2}}{2}\right)$.

###  4.3 Kernels for categorical variables

Let us study a kernel
$k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{V}\subset\mathbb{R}$, where
$\mathcal{X}$ is a finite set of $C$ possible values (categories). We can
encode these categories numerically as integers
$\mathcal{X}=\\{1,\ldots,C\\}$. Because there are only $C^{2}$ possible input
combinations $(v,w)$ for $k$, and therefore $\lvert\mathbb{V}\rvert\leq
C^{2}$, we can list them in the $C\times C$ matrix $\bm{C}$ which has elements
$\left[\bm{C}\right]_{v,w}=k(v,w)$. If $k$ is symmetric, the symmetric square
matrix $\bm{C}$ has the orthogonal eigendecomposition

$\bm{C}=\bm{\Theta}\bm{D}\bm{\Theta}^{\top},$ (16)

where $\bm{D}$ is the diagonal matrix containing the eigenvalues $d_{c}$,
$c=1,\ldots,C$ on the diagonal and $\bm{\Theta}$ has the $C$ eigenvectors as
its columns. For each column $c$, we can define function
$\varphi_{c}:\mathcal{X}\rightarrow\mathbb{R}$ so that
$\varphi_{c}(v)=\left[\bm{\Theta}\right]_{v,c}$. We see that

$\displaystyle
k(v,w)=\left[\bm{C}\right]_{v,w}=\left[\bm{\Theta}\bm{D}\bm{\Theta}^{\top}\right]_{v,w}$
(17)
$\displaystyle=\sum_{c=1}^{C}d_{c}\left[\bm{\Theta}\right]_{v,c}\left[\bm{\Theta}\right]_{w,c}=\sum_{c=1}^{C}d_{c}\varphi_{c}(v)\varphi_{c}(w),$
(18)

meaning that we have written $k$ in the form of Eq. 11 with $M=C$ and
$\psi_{m}(\bm{x})=\sqrt{d_{m}}\varphi_{m}(\bm{x})$. Note that this is an exact
function decomposition for $k$ and not an approximation. The complexity of
computing the eigendecomposition is $\mathcal{O}(C^{3})$, but in typical
applications $C\ll N$ and this is not a bottleneck. Actually, for example for
the ZS kernel and other CS kernels, the eigenvalues have a closed form and the
corresponding eigenbasis is known (see Suppl. Section 2). Furthermore, if $k$
does not depend on any hyperparameters, the eigendecomposition only needs to
be done once before parameter inference. If it is of type
$k(\bm{x},\bm{x}^{\prime})=\alpha^{2}k_{0}(\bm{x},\bm{x}^{\prime})$ where
$\alpha$ is the only parameter, the decomposition can obviously be done just
for $k_{0}$ which again has no parameters. Evaluating functions $\varphi_{m}$
is easy as it corresponds to just looking up a value from the matrix
$\bm{\Theta}$.

###  4.4 Mixed-domain product kernels

We now consider approximating a product kernel
$k(\bm{x},\bm{x}^{\prime})=\prod_{i=1}^{P}k_{i}(\bm{x},\bm{x}’)$ with
$\tilde{k}(\bm{x},\bm{x}^{\prime})=\prod_{i=1}^{P}\tilde{k}_{i}(\bm{x},\bm{x}’)$,
where for each $\tilde{k}_{i}$ we have an available decomposition

$\tilde{k}_{i}(\bm{x},\bm{x}^{\prime})=\sum_{m=1}^{M_{i}}\psi_{i,m}(\bm{x})\psi_{i,m}(\bm{x}^{\prime}),$
(19)

which might be an approximation or an exact decomposition of $k_{i}$. The
total approximation is

$\displaystyle\tilde{k}(\bm{x},\bm{x}^{\prime})=\prod_{i=1}^{P}\sum_{m=1}^{M_{i}}\psi_{i,m}(\bm{x})\psi_{i,m}(\bm{x}^{\prime})$
(20)
$\displaystyle=\sum_{m_{1}=1}^{M_{1}}\ldots\sum_{m_{P}=1}^{M_{P}}\psi_{m_{1},\ldots,m_{P}}^{*}(\bm{x})\psi_{m_{1},\ldots,m_{P}}^{*}(\bm{x}^{\prime})$
(21)

where
$\psi^{*}_{m_{1},\ldots,m_{P}}(\bm{x})=\prod_{i=1}^{P}\psi_{i,m_{i}}(\bm{x})$.
We have now a representation of the product kernel in the form of Eq. 11 with
$M=\prod_{i=1}^{P}M_{i}$ sum terms. Note that since the individual kernels in
the product kernel can be both categorical and continuous, Eq. 21 provides a
kernel representation for mixed-domain GPs with product kernels. Also note
that $M$ grows exponentially with $P$.

###  4.5 Mixed-domain sum kernels

The most general kernels that we consider are of the form
$k(\bm{x},\bm{x}^{\prime})=\sum_{j=1}^{J}\prod_{i=1}^{P_{j}}k_{i,j}(\bm{x},\bm{x}’)$,
where $P_{j}\geq 1$ is the number of product factors in the sum term $j$. If
each $k_{i,j}$ has a (possibly approximate) decomposition

$\tilde{k}_{i,j}(\bm{x},\bm{x}^{\prime})=\sum_{m=1}^{M(i,j)}\psi_{i,j,m}(\bm{x})\psi_{i,j,m}(\bm{x}^{\prime}),$
(22)

with $M(i,j)$ sum terms, we can approximate $k(\bm{x},\bm{x}^{\prime})$ with

$\displaystyle\tilde{k}(\bm{x},\bm{x}^{\prime})=\sum_{j=1}^{J}\prod_{i=1}^{P_{j}}\tilde{k}_{i,j}(\bm{x},\bm{x}’)$
(23)
$\displaystyle=\sum_{j=1}^{J}\sum_{m_{1}=1}^{M(P_{1},1)}\ldots\sum_{m_{P_{j}}=1}^{M(P_{j},j)}\psi_{j,m_{1},\ldots,m_{P_{j}}}^{*}(\bm{x})\psi_{j,m_{1},\ldots,m_{P_{j}}}^{*}(\bm{x}^{\prime})$

where
$\psi^{*}_{j,m_{1},\ldots,m_{P_{j}}}(\bm{x})=\prod_{i=1}^{P_{j}}\psi_{i,j,m_{i}}(\bm{x})$
for each $j$. Now we have a sum representation (Eq. 11) of the kernel with
$M=\sum_{j=1}^{J}\prod_{i=1}^{P_{j}}M(i,j)$ terms.

###  4.6 Mixed kernels for longitudinal data

In our framework, we consider mixed kernels
$k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$, where $\mathcal{X}$ is
a mixed space of both continuous and categorical dimensions, consisting of
multiplication and addition so that

$k(\bm{x},\bm{x}^{\prime})=\sum_{j=1}^{J}\left(\prod_{q=1}^{Q_{j}}k_{j,q}(\bm{x},\bm{x}^{\prime})\prod_{r=1}^{R_{j}}l_{j,r}(\bm{x},\bm{x}^{\prime})\right),$
(24)

where each $k_{j,q}:\mathbb{R}\times\mathbb{R}$ is isotropic and depends only
one continuous dimension of $\bm{x}$ and each
$l_{j,r}:\\{1,\ldots,C_{j,r}\\}\rightarrow\mathbb{R}$ depends only on one
categorical dimension of $\bm{x}$, which has $C_{j,r}$ different categories.
For each $k_{j,q}$, we use the basis function approximation (Eq. 13) with
$B_{j,q}$ basis functions and domain $\Omega=[L_{j,q},L_{j,q}]$, and for each
$l_{j,r}$ the exact decomposition (Eq. 18). Using Eq. 23, we can write
$k(\bm{x},\bm{x}^{\prime})$ in the format of Eq. 11 with

$M=\sum_{j=1}^{J}\left(\prod_{q=1}^{Q_{j}}B_{j,q}\prod_{r=1}^{R_{j}}C_{j,r}\right)$
(25)

terms. In each term, the function $\psi(\bm{x})$ is a product of $Q_{j}$
factors $\sqrt{s}\phi(\bm{x})$ and $R_{j}$ factors $\sqrt{d}\varphi(\bm{x})$.

As an example, if $Q_{j}=R_{j}=1$ for each $j$ and we use $B_{j,1}=B$ basis
functions for all components, then the scalability is $\mathcal{O}(NM+M)$
where $M=B\sum_{j=1}^{J}C_{j,1}$. Further, if each categorical variable has
$C$ many different values, then the scalability is $\mathcal{O}(NM+M)$, where
$M=BJC$.

##  5 RESULTS

We demonstrate the scalability and accuracy of the presented approach using
experiments with simulated and real data. In all experiments we use the
Dynamic HMC algorithm of Stan (version 2.27) (Carpenter et al., 2017), with
target acceptance rate set to $0.95$, for MCMC sampling the parameters of our
approximate models111Code will be made available at
https://github.com/jtimonen/scalable-mixed-domain-GPs.. All models are fitted
by running four independent MCMC chains for 2000 iterations each, discarding
the first half of each chain as warmup. In all experiments, we use Student-$t$
priors with $20$ degrees of freedom for the $\alpha$ kernel parameters, log-
normal priors with mean 0 and scale 1 for kernel lengthscale parameters
$\ell$. In Experiments 1 and 2 the noise variance parameter $\sigma^{2}$ of
the Gaussian observation model has an Inverse-Gamma prior with shape
parameters $\alpha=1$ and $\beta=2$. Priors are on normalized data scale,
meaning that continuous variables are standardized to zero mean and unit
variance during inference.

In all experiments we use the same number of basis functions $B$ for each
approximate continuous kernel. We use also for all approximate components the
same domain scaling factor $c$, which is defined so that $\Omega=[L,L]$ with
$L$ being $c$ times the half-range of the continuous covariate of the
approximated kernel (Riutort-Mayol et al., 2020).

Experiments 1-2 are run on a modern CentOS 7 computing cluster and Experiment
3 on a laptop computer.

Figure 1: The posterior predictive mean for four different approximate models
in one replication of Experiment 1. The yellow line shows the posterior
predictive mean for the corresponding exact GP model as a reference. Black
dots are training data and red crosses test data.

###  5.1 Experiment 1: Simulation study

In the first experiment we create simulated longitudinal data consisting of
categorical variables $id$ and $z$, and a continuous variable $age$. We create
data with 9 individuals, where individuals with $id\in\\{1,4,7\\}$ belong to
group $z=1$, individuals with $id\in\\{2,5,8\\}$ belong to group $z=2$ and
individuals $id\in\\{3,6,9\\}$ belong to group $z=3$. For each individual 1-6,
we create $\frac{N_{\text{train}}}{6}$ observations at time points where $age$
is drawn uniformly from the interval $[0,10]$, and $N_{\text{train}}$ is
varied as $\\{60,90,120,150,180,210\\}$. For individuals 7-9,
$\frac{N_{\text{test}}}{3}$ observations are created similarly, with
$N_{\text{test}}=150$.

Figure 2: Mean log predictive density for test data in Experiment 1. Black
dashed line corresponds to the exact model. Results are averages over 30
replications of the experiment.

We consider an additive GP model $f=f^{(1)}(age)+f^{(2)}(age,z)$ with kernels

$\begin{cases}k_{1}&=\alpha_{1}^{2}k_{\text{EQ}}(age,age^{\prime}\mid\ell_{1})\\\
k_{2}&=\alpha_{2}^{2}k_{\text{ZS}}(z,z)k_{\text{EQ}}(age,age^{\prime}\mid\ell_{2})\end{cases}$
(26)

and we simulate a realization of
$\bm{f}\in\mathbb{R}^{N_{\text{test}}+N_{\text{train}}}$ data using
$\alpha_{1}=\alpha_{2}=1$, $\ell_{1}=2$ and $\ell_{2}=1$. We then generate a
response variable measurements $\bm{y}=100+10\cdot(\bm{f}+\bm{\epsilon})$,
where $\epsilon_{i}\sim\mathcal{N}(0,0.5^{2})$ and the realization $\bm{f}$
represents the ground truth signal.

Data from individuals 1-6 is used in training, while data from individuals 7-9
is left for testing. Using the training data, we fit an exact and approximate
model with the correct covariance structure from Eq. 26 using the Gaussian
likelihood model. Exact model is fitted with lgpr (Timonen et al., 2021),
which also uses Stan for MCMC. The exact model utilizes the marginalization
approach for GPs, as Gaussian observation model is specified.

Figure 1 shows the posterior predictive mean of the exact model and different
approximate models with $c=1.5$, using $N_{\text{train}}=150$. We see that
with $B=16$ and $B=32$ the mean predictions are indistinguishable from the
exact model. We fit the approximate model using different values of $B$ and
$c$, and repeat the experiment using different values for $N_{\text{train}}$.
Results in Figure 3 validate empirically that the runtime scales linearly as a
function of both $N_{\text{train}}$ and $B$.

We compute the mean log predicive density (MLPD), at test points
$\bm{y}^{*}\in\mathbb{R}^{150}$ (see Suppl. Section 3 for details about out-
of-sample prediction and MLPD). Results in Figure 2 show that the MLPD of the
approximate model approaches the exact model as $B$ grows. It is seen that
with small data sizes and small $B$, the predictive performance can actually
be better than that of the exact model, possibly because the coarser
approximation is a simpler model that generalizes better in this case.

Figure 3: Runtimes of fitting the exact and approximate models in Experiment
1. a) Exact model vs. approximation with $B=64$. b) Approximations using
different values for $B$. The markers show the average time taken to run a
chain, when a total of $30\times 4$ MCMC chains were run for 2000 iterations
each. The vertical error bars show $\pm$ one standard deviation (not shown for
approximate model in panel a). Note the smaller y-axis scale in panel b. We
see empirically that the runtime of the approximate model scales linearly in
both $N$ and $B$.

###  5.2 Experiment 2: Canadian weather data

We analyse data that consists of yearly average temperature measurements in 35
Canadian weather stations (Ramsay and Silverman, 2005). There is a total of
$N=35\times 365=12775$ data points, which are daily temperatures at the 35
locations, averaged over the years 1960-1994. We fit an additive GP model
$f=f^{(1)}(day)+f^{(2)}(day,region)+f^{(3)}(day,station)$, with Gaussian
likelihood using the EQ kernel for $f^{(1)}$ and the product EQ$\times$ZS
kernel for $f^{(2)}$ and $f^{(3)}$.

Figure 4: Results for the Canadian weather data experiment with $B=32$ and
$c=1.5$. Figures a)-c) show the marginal posterior distribution of each of the
three components $f^{(1)}$, $f^{(2)}$ and $f^{(3)}$ (mean $\pm$ two times
standard deviation). Standard deviation is not shown for $f^{(3)}$ for
clarity. The functions are on the standardized scale (response variable
normalized to zero mean and unit variance). We see for example that the
regions tend to have larger differences during winter.

We used domain scaling factor $c=1.5$ for all components and ran the 4 MCMC
chains in parallel using 4 CPU cores. This was repeated with different values
of $B=8,12,16,24,32$, where $B$ is the number of basis functions for each
components. Total runtimes for fitting the models were $4.10,7.24,7.47,14.18$
and $18.76$ hours, respectively. The posterior distributions of each model
component with $B=32$ are in Figure 4. The posterior predictive distribution
for each station separately is visualized in Suppl. Figure 1.

###  5.3 Experiment 3: US presidential election prediction

In the last example we demonstrate a beta-binomial observation model and model
the vote share of the Republican Party in each state in US presidential
elections. By two-party vote share we mean proportion of votes cast to the
Republican candidate divided by the sum of votes to both the Republican and
Democratic candidates222Data is from (MIT Election Data and Science Lab,
2017).. Following Trangucci (2017), Washington DC is excluded from the
analysis. We use data from the 1976-2016 elections as training data, meaning
that $N=50\times 11=550$.

We fit an additive GP model
$f=f^{(1)}(year)+f^{(2)}(year,region)+f^{(3)}(year,state)$, with beta-binomial
observation model using the EQ kernel for $f^{(1)}$ and the product
EQ$\times$ZS kernel for $f^{(2)}$ and $f^{(3)}$. The observation model is

$\text{y}^{R}\sim\text{Beta-Binomial}(\text{y}^{R}+\text{y}^{D},a,b),$ (27)

where $y^{R},y^{D}$ are the number of votes for the Republican and Democratic
parties, respectively, $a=\rho(\gamma^{-1}-1)$, $b=(1-\rho)(\gamma^{-1}-1)$,
and $\rho=\text{inv-logit}(f+w_{0})$. We use a $\mathcal{N}(0,0.5^{2})$ prior
for the intercept $w_{0}$ and a Log-Normal(1,1) prior for the $\gamma$
parameter.

Fitting a model with $B=24$ and $c=2.0$ on a 2018 MacBook Pro computer (2.3
GHz Quad-Core Intel i5 CPU), running the 4 chains in parallel, took
approximately 18 minutes. The posterior distributions of each model component
are in Figure 5. See also Supplementary Figure 2, where we have also
visualized the data from 2020 election to validate that the model predicts
well into the future.

##  6 CONCLUSION

Gaussian processes offer an attractive framework for specifying flexible
models using a kernel language. The computational cost of their exact
inference however limits possible applications to small data sets. Our
scalable framework opens up a rich class of GP models to be used in large
scale applications of various fields of science as the computational
complexity is linear with respect to data size. We have presented a scalable
approximation scheme for mixed-domain covariance functions, and demonstrated
its use in the context of Bayesian GP regression. However, it can also be
applied in GP applications where the kernel hyperparameters are optimized
using a marginal likelihood criterion.

We recall that we have assumed that the categorical kernels are symmetric, and
continuous kernels are stationary. Non-stationary effects can still be modeled
by applying a warping on the input first, and then using a stationary kernel
(see for example (Cheng et al., 2019)). Another limitation of the approach is
that when the number of product terms in a kernel grows, the total number of
basis functions required for that component grows exponentially and can become
too large. This still leaves us with a large class of mixed-domain GP models
that are scalable.

Figure 5: Results for the US election prediction experiment with $B=24$ and
$c=2.0$. Figures a)-c) show the marginal posterior distribution of each of the
three components $f^{(1)}$, $f^{(2)}$ and $f^{(3)}$ (mean $\pm$ two times
standard deviation). Standard deviation is not shown for $f^{(2)}$ and
$f^{(3)}$ for clarity.

#### Acknowledgements

We thank Aki Vehtari and Gleb Tikhonov for useful comments on draft versions
of this manuscript, and acknowledge the computational resources provided by
Aalto Science-IT, Finland. This work was supported by the Academy of Finland
and Bayer Oy.

## References

  * Cao et al. (2015) Y. Cao, M. A. Brubaker, D. J. Fleet, and A. Hertzmann.  Efficient optimization for sparse Gaussian process regression.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 37(12):2415–2427, 2015.
  * Carpenter et al. (2017) B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell.  Stan: A probabilistic programming language.  _Journal of Statistical Software_ , 76(1):1–32, 2017.
  * Cheng et al. (2019) L. Cheng, S. Ramchandran, T. Vatanen, N. Lietzen, R. Lahesmaa, A. Vehtari, and H. Lähdesmäki.  An additive Gaussian process regression model for interpretable non-parametric analysis of longitudinal data.  _Nature Communications_ , 10, 2019.
  * Chung et al. (2020) I. Chung, S. Kim, J. Lee, K. J. Kim, S. J. Hwang, and E. Yang.  Deep mixed effect model using Gaussian processes: A personalized and reliable prediction for healthcare.  _The Thirty-Fourth AAAI Conference on Artificial Intelligence (AAAI-20)_ , 2020.
  * Deng et al. (2017) X. Deng, C. D. Lin, K.-W. Liu, and R. K. Rowe.  Additive Gaussian process for computer models with qualitative and quantitative factors.  _Technometrics_ , 59(3):283–292, 2017.
  * Fortuin et al. (2021) V. Fortuin, G. Dresdner, H. Strathmann, and G. Rätsch.  Sparse Gaussian processes on discrete domains.  _IEEE Access_ , 9:76750–76758, 2021.
  * Garrido-Merchán and Hernández-Lobato (2020) E. C. Garrido-Merchán and D. Hernández-Lobato.  Dealing with categorical and integer-valued variables in Bayesian optimization with Gaussian processes.  _Neurocomputing_ , 380:20–35, 2020.
  * Kaufman and Sain (2010) C. G. Kaufman and S. R. Sain.  Bayesian functional ANOVA modeling using Gaussian process prior distributions.  _Bayesian Analysis_ , 5(1):123–149, 2010.
  * Kimeldorf and Wahba (1970) G. S. Kimeldorf and G. Wahba.  A Correspondence Between Bayesian Estimation on Stochastic Processes and Smoothing by Splines.  _The Annals of Mathematical Statistics_ , 41(2):495 – 502, 1970.
  * Liu et al. (2020) H. Liu, Y.-S. Ong, X. Shen, and J. Cai.  When Gaussian process meets big data: A review of scalable GPs.  _IEEE Transactions on Neural Networks and Learning Systems_ , 31(11):4405–4423, 2020.
  * MIT Election Data and Science Lab (2017) MIT Election Data and Science Lab.  U.S. President 1976–2020, V6, 2017.
  * Pedersen et al. (2019) E. Pedersen, D. Miller, G. Simpson, and N. Ross.  Hierarchical generalized additive models in ecology: An introduction with mgcv.  _PeerJ_ , (5), 2019.
  * Quiñonero Candela and Rasmussen (2005) J. Quiñonero Candela and C. E. Rasmussen.  A unifying view of sparse approximate Gaussian process regression.  _J. Mach. Learn. Res._ , 6:1939–1959, 2005.
  * Quintana et al. (2016) F. A. Quintana, W. O. Johnson, L. E. Waetjen, and E. B. Gold.  Bayesian nonparametric longitudinal data analysis.  _Journal of the American Statistical Association_ , 111(515):1168–1181, 2016.
  * Ramsay and Silverman (2005) J. Ramsay and B. W. Silverman.  _Functional Data Analysis_.  Springer, New York, NY, 2nd edition, 2005.
  * Rasmussen and Williams (2006) C. E. Rasmussen and C. K. I. Williams.  _Gaussian Processes for Machine Learning_.  MIT Press, Cambridge, Massachusetts, 2006.
  * Riutort-Mayol et al. (2020) G. Riutort-Mayol, P.-C. Bürkner, M. R. Andersen, A. Solin, and A. Vehtari.  Practical hilbert space approximate Bayesian Gaussian processes for probabilistic programming, 2020.
  * Roustant et al. (2020) O. Roustant, E. Padonou, Y. Deville, A. Clément, G. Perrin, J. Giorla, and H. Wynn.  Group kernels for Gaussian process metamodels with categorical inputs.  _SIAM/ASA Journal on Uncertainty Quantification_ , 8(2):775–806, 2020.
  * Snelson and Ghahramani (2006) E. Snelson and Z. Ghahramani.  Sparse Gaussian processes using pseudo-inputs.  In Y. Weiss, B. Schölkopf, and J. Platt, editors, _Advances in Neural Information Processing Systems_ , volume 18. MIT Press, 2006.
  * Solin and Särkkä (2019) A. Solin and S. Särkkä.  Hilbert space methods for reduced-rank Gaussian process regression.  _Statistics and Computing_ , 2019.
  * Timonen et al. (2021) J. Timonen, H. Mannerström, A. Vehtari, and H. Lähdesmäki.  lgpr: an interpretable non-parametric method for inferring covariate effects from longitudinal data.  _Bioinformatics_ , 37(13):1860–1867, 2021.
  * Titsias (2009) M. Titsias.  Variational learning of inducing variables in sparse Gaussian processes.  In D. van Dyk and M. Welling, editors, _Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics_ , volume 5 of _Proceedings of Machine Learning Research_ , pages 567–574, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 2009. PMLR.
  * Trangucci (2017) R. Trangucci.  Hierarchical Gaussian processes in Stan, 2017.  URL https://mc-stan.org/events/stancon2017-notebooks/stancon2017-trangucci-hierarchical-gps.pdf.
  * Verbeke and Molenberghs (2000) G. Verbeke and G. Molenberghs.  _Linear Mixed Models for Longitudinal Data_.  Springer, New York, NY, 2000.
  * Wang et al. (2021) L. Wang, S. Yerramilli, A. Iyer, D. Apley, P. Zhu, and W. Chen.  Scalable Gaussian processes for data-driven design using big data with categorical factors.  _arXiv:2106.15356_ , 2021.
  * Wilson and Nickisch (2015) A. G. Wilson and H. Nickisch.  Kernel interpolation for scalable structured Gaussian processes (KISS-GP).  _CoRR_ , abs/1503.01057, 2015.  URL http://arxiv.org/abs/1503.01057.
  * Zhang and Notz (2015) Y. Zhang and W. I. Notz.  Computer experiments with qualitative and quantitative variables: A review and reexamination.  _Quality Engineering_ , 27:2–13, 2015.
  * Zhang et al. (2020) Y. Zhang, S. Tao, W. Chen, and D. W. Apley.  A latent variable approach to gaussian process modeling with qualitative and quantitative factors.  _Technometrics_ , 62(3):291–302, 2020.


n}{n}}\right)\frac{n}{k}\leqslant\left\lVert
c_{p}(\mathbf{z})\right\rVert^{2}\leqslant\left(1+\sqrt{\frac{400k\log
n}{n}}\right)\frac{n}{k}\,,$ hence the property follows by a union bound. ∎

Next we prove 4.5.

###### Proof.

For each $\ell\in[t]\,,$ define

$\operatorname*{\mathbb{E}}\left[\hat{\mathbf{X}}(\mathbf{G}_{\ell})_{ij}\;\middle|\;\mathbf{f}_{\ell}(\mathbf{z})_{i}=\mathbf{f}_{\ell}(\mathbf{z})_{j}\right]=:C^{*}_{\ell}$

and let

$C^{*}=\sum_{\ell\in[t]}C_{\ell}\,.$

By independence of the observations, the inequality of 4.5 for the case case
$\mathbf{z}_{i}=\mathbf{z}_{j}$ follows with an application of Hoeffding’s
inequality.

The case $\mathbf{z}_{i}\neq\mathbf{z}_{j}$ needs a bit more work. By Theorem
3.1, for each $\ell\in[t]$, we have

$\operatorname*{\mathbb{E}}\left[\hat{\mathbf{X}}(\mathbf{G}_{\ell})_{ij}\mathrel{}\middle|\mathrel{}\mathbf{f}_{\ell}(\mathbf{z})_{i}=\mathbf{f}_{\ell}(\mathbf{z})_{j}\right]=C^{*}_{\ell}\,,$

and

$\operatorname*{\mathbb{E}}\left[\hat{\mathbf{X}}(\mathbf{G}_{\ell})_{ij}\mathrel{}\middle|\mathrel{}\mathbf{f}_{\ell}(\mathbf{z})_{i}\neq\mathbf{f}_{\ell}(\mathbf{z})_{j}\right]\leqslant
C^{*}_{\ell}-C_{\ell}\,.$

Therefore,

$\displaystyle\operatorname*{\mathbb{E}}\left[\hat{\mathbf{X}}(\mathbf{G}_{\ell})_{ij}\mathrel{}\middle|\mathrel{}\mathbf{z}_{i}\neq\mathbf{z}_{j}\right]=$
$\displaystyle\operatorname*{\mathbb{E}}\left[\hat{\mathbf{X}}(\mathbf{G}_{\ell})_{ij}\mathrel{}\middle|\mathrel{}\mathbf{f}_{\ell}(\mathbf{z})_{i}=\mathbf{f}_{\ell}(\mathbf{z})_{j}\right]\cdot\mathbb{P}\left(\mathbf{f}_{\ell}(\mathbf{z})_{i}=\mathbf{f}_{\ell}(\mathbf{z})_{j}\mathrel{}\middle|\mathrel{}\mathbf{z}_{i}\neq\mathbf{z}_{j}\right)$
$\displaystyle+\operatorname*{\mathbb{E}}\left[\hat{\mathbf{X}}(\mathbf{G}_{\ell})_{ij}\mathrel{}\middle|\mathrel{}\mathbf{f}_{\ell}(\mathbf{z})_{i}\neq\mathbf{f}_{\ell}(\mathbf{z})_{j}\right]\cdot\mathbb{P}\left(\mathbf{f}_{\ell}(\mathbf{z})_{i}\neq\mathbf{f}_{\ell}(\mathbf{z})_{j}\mathrel{}\middle|\mathrel{}\mathbf{z}_{i}\neq\mathbf{z}_{j}\right)$
$\displaystyle\leqslant$ $\displaystyle
C^{*}_{\ell}-C_{\ell}\cdot\frac{1}{2}\,,$

and so

$\operatorname*{\mathbb{E}}\left[\sum_{\ell\in[t]}\hat{\mathbf{X}}(\mathbf{G}_{\ell})_{ij}\;\middle|\;\mathbf{z}_{i}\neq\mathbf{z}_{j}\right]\leqslant\sum_{\ell\in[t]}\left(C^{*}_{\ell}-C_{\ell}\cdot\frac{1}{2}\right)=C^{*}-\frac{\overline{C}\cdot
t}{2}\,.$

The inequality of 4.5 for case $\mathbf{z}_{i}\neq\mathbf{z}_{j}$ follows with
another application of Hoeffding’s inequality. ∎

Now we prove Lemma 4.8.

###### Proof of Lemma 4.8.

Let $A_{i}$ be a $(p,q)$-representative and $A_{j}$ be a
$(p^{\prime},q)$-representative, for $p\,,p^{\prime}\in[k]\,.$ It suffices to
show that if $p=p^{\prime}$ then $\left\lVert
A_{i}-A_{j}\right\rVert^{2}\leqslant n/k$ and otherwise $\left\lVert
A_{i}-A_{j}\right\rVert^{2}>n/k\,.$ Then no $q$-representative index remains
unassigned at the end of step $1$. By the reverse triangle inequality

$\displaystyle\left\lvert\left\lVert A_{i}-A_{j}\right\rVert-\left\lVert
c_{p}(z)-c_{p^{\prime}}(z)\right\rVert\right\rvert$
$\displaystyle\leqslant\left\lVert
A_{i}-A_{j}-c_{p}(z)+c_{p^{\prime}}(z)\right\rVert$
$\displaystyle\leqslant\left\lVert A_{i}-c_{p}(z)\right\rVert+\left\lVert
A_{j}-c_{p^{\prime}}(z)\right\rVert$ $\displaystyle\leqslant 2n\cdot
e^{-q\cdot\overline{C}^{2}\cdot t}\,.$

For $p=p^{\prime}$ we have $\left\lVert
c_{p}(z)-c_{p^{\prime}}(z)\right\rVert=0$ and by choice of $t\,,$ the first
inequality follows. For $p\neq p^{\prime}$ we have $\left\lVert
c_{p}(z)-c_{p^{\prime}}(z)\right\rVert\geqslant\frac{n}{k}(2-o(1))$ since $z$
is balanced and the second inequality follows again by choice of $t$. ∎