\section{Introduction and Related Work}
Most of the recent interest in the area of \textbf{word embeddings} was triggered by the Word2Vec algorithm proposed in \cite{Mikolov2013}, which provided an efficient way to learn word embeddings by predicting words based on their context and using negative sampling. Word embeddings have become the usual input to natural language processing (NLP) tasks, but also tasks for which previously knowledge graphs were being used. Applications range from text classification \cite{Kim2014} to machine translation \cite{Kalchbrenner2013,Sutskever,Cho2014,Bahdanau2014}, question answering \cite{Khot,Seo,Parikh} and knowledge graph construction and completion\cite{fu2014semHierarchiesViaEmbeddings}. However, despite recent efforts \cite{Li,Garcia2018,Zeiler2014}, the nature and extent of the knowledge captured by such embeddings and how they contribute to accomplish the goal in question is still hard to interpret.

Embeddings have shown the ability to learn relations between words. However, most benchmarks \cite{Schnabel2015} focus on a specific family of relations involving relations similarity and analogy. It is unclear whether such relations simply happen to be well aligned with the statistic analysis involved in the computation of the embeddings, or whether the embeddings are capable of capturing a wider range of relational knowledge. Indeed, standard benchmarks show evidence that word embeddings may capture more types of relations. It has been shown that algorithms like FastText \cite{Bojanowski2017}, GloVe \cite{pennington2014glove} and Swivel \cite{Shazeer2016} learn embeddings that capture lexical and semantic information. However, the current lack of a standard evaluation practice make it hard to study which specific relations embeddings can effectively capture, nor how to best quantify the signal of such relations.

At the same time, embeddings as a knowledge representation mechanism is being explored by the traditionally symbolic community that produced semantic networks and knowledge graphs. Algorithms based on \textbf{knowledge graph (KG) embeddings}, like RDF2Vec\cite{Ristoski2016}, ProjE\cite{Shi2016}, TransE \cite{Bordes2013} and HolE\cite{nickel2016HolE}, learn embeddings representing the concepts, words and relations contained in a KG and provide a vector representation of the knowledge that is explicitly described in it. However, such KG embeddings may only encode knowledge that is already represented in the KG. One application of embeddings for this community is \textbf{generic entity-based KG completion and refinement} \cite{riedel2013univSchemas,melo2017relErrorsInKGs,nickel2016relMLForKGs,Paulheim2017KGRefinement,lin2016relExtractOverInstances}, which tries to use large text corpora to complete a partial knowledge graph. However, such efforts have focused on encyclopedic (e.g. DBpedia) or domain-specific (family, commerce, finance, law) relations between entities rather than lexical relations. On the opposite direction, KGs are also used to refine vector space representations as in \cite{Faruqui2015}. Finally, efforts trying to learn \textbf{lexical semantic relations}~\cite{shwartz2016LexNet,gabor2017vecs4semrels,turney2015recognizingLexEntailment,roller2014distribHypernymDetection,fu2014semHierarchiesViaEmbeddings} have started looking at whether word embeddings (and similar distributional approaches) can be used to predict certain types of lexical semantic relations. Many of these approaches suffer from the difficulty to generate a training dataset that serves this purpose~\cite{levy2015doEmbedsLearnRels} and most of them focus on hypernymy relations or lexical inference, still a limited fragment of the whole spectrum of possible lexical semantic relations. Also, WordNet tends to be the only source of evidence used in such work, which may hinder reaching a general understanding of the matter.

\begin{figure*}[h]
  \includegraphics[width=0.9\textwidth]{rel-prediction-flat}
  \caption{Generic approach for evaluating the predictive power of word and concept embeddings using a knowledge graph as a silver standard.}
 \label{fig:relpredict}
\end{figure*}

In previous work\cite{denaux2017vecsigrafo,denauxvecsigrafo}\footnote{\href{https://bit.ly/2LykQ9T}{Vecsigrafo: Corpus-based Word-Concept Embeddings}} on joint word-concept embeddings, these outperformed word-only and knowledge graph approaches over a selection of 14 benchmarks on semantic similarity and relatedness as well as word-concept and hypernym relation prediction tasks. In doing so, we suggested it may be possible to extract knowledge about lexical semantic relations from (joint word-concept) embeddings. In this paper, we take into account lessons from the various communities discussed above to propose a method for evaluating to what extent pre-trained embeddings contain knowledge about a wide ranging of relations encoded in a KG. Our contributions are: (i) describing a (largely automated) method for (i.a) word/concept relational dataset generation from a KG (i.b) training of machine learning models and (i.c) evaluation of the generated datasets and trained models; (ii) application of this method to analyse to what extent (ii.a) various corpus- and KG-based pre-trained embeddings capture (ii.b) various types of lexico-semantic relational knowledge and (ii.b) what the effect is various factors such as the size of the corpus and the type of dataset (word or concept).

In Section~\ref{sec:genericMethod} we describe a generic method for measuring the predictive power of embeddings for specific relations. In section~\ref{sec:lexicoSemanticStudy} we describe how we apply this methodology to study lexico-semantic relations and in section~\ref{sec:analysis} we analyse the gathered data.

\section{Measuring Relation Predictive Power of Embeddings with a KG}
\label{sec:genericMethod}

In this section we propose a method for studying the predictive power of word and concept embeddings, depicted in Figure~\ref{fig:relpredict}. The goal is to study how much knowledge about relations is encoded in embedding spaces and can be exploited by machine learning models. Note that our goal is \emph{not} to generate the best models for predicting relations (which could be achieved by combining different sources of evidence). The method consists of three main phases: dataset generation, model training and prediction results analysis.

\paragraph{Preliminaries} We define a \textbf{Knowledge Graph} as a tuple $\langle N, R, E \rangle$, where $N = C \cup I \cup W$ is a set of node identifiers, typically referring to concepts $c \in C$, instances $i \in I$ or words $w \in W$ (human readable names for concepts and instances); $R$ is a set of relation types and $E$ is a set of triples of the form $(r\ n_i\ n_j)$ where $r \in R$ and $n_i,n_j \in N$. We define an \textbf{embedding space} $S$ as a tuple $\langle V, F_d \rangle$, where $V$ is a vocabulary, and $F_d: V \longrightarrow \mathbb{R}^d$ is a function that maps elements in the vocabulary to its embedding: a vector of dimension $d$.

The main \emph{inputs} for our method are a KG $k = \langle N_k, R_k, E_k \rangle$ and one or more embedding spaces $s_i \in S'_k \subset S_k \subset S$. $S_k$ is the set of all embedding spaces that have a vocabulary that overlaps with $N_k$. $S'_k$ is a subset of embedding spaces to study. The main \emph{output} of our method is a classification of tuples $(r\ s)$, where $r \in R_k$ and $s \in S'_k$, into \emph{predictable} and \emph{non-predictable}. Furthermore, for each such tuple, we also derive absolute and relative relation prediction metrics, providing a numerical assessment of how well the embedding space $s$ captures or encodes knowledge about relationship $r$.

\paragraph{Dataset generation}
In the first phase, we aim to generate datasets $\delta \in D$, where each dataset $\delta$ is a finite set of triples of the form $\delta_r = \{\langle n_i\ n_j\ l\rangle\}$, where $l$ is a classification label such that $l = 1$ if $(r\ n_i\ n_j) \in E_k$ and $l = 0$ otherwise. We generate datasets between \emph{words} $\delta^w$ with tuples $\langle w_i\ w_j\ l\rangle$, word-concepts $\delta^{wc}$ with tuples $\langle w_i\ c_j\ l\rangle$ and concepts $\delta^c$ with tuples $\langle c_i\ c_j\ l\rangle$.

As a first step, we select a \textbf{seed vocabulary} $V_{\mathrm{seed}} = \bigcap_{i=1}^{|S'_k|} V_{s_i}$, as the intersection of all the studied vocabularies.
Next, $V_{\mathrm{seed}}$ and $k$ are used to \textbf{extract a set of positive pairs} for each $r \in R_k$; these correspond to partial datasets $\delta_r^+ = \{\langle n_i\ n_j\  1 \rangle \mid (r\ n_i\ n_j) \in E_k \wedge n_i,n_j \in V_{\emph{seed}}\}$. KGs frequently define relations at a concept level, therefore, this initial dataset will typically only contain concepts. To also generate word-concept datasets for such relations, we also extract partial datasets
 $\{ \langle w_i\ c_j\ 1\rangle \mid (r\ c_i\ c_j) \wedge (r_w\ w_i\ c_i) \}$, where $r_w$ is the word-to-concept relation in $k$. Similarly, we can extract a word dataset for $r$ using  $\{ \langle w_i\ w_j\ 1\rangle \mid (r\ c_i\ c_j) \wedge (r_w\ w_i\ c_i) \wedge (r_w\ w_j\ c_j) \}$.

  Besides generating positive datasets for relations $r \in R_k$, we also generate \emph{random datasets}. These datasets serve as baselines and will be used later on to prune biased datasets. Since the generated datasets $\delta_r^+$ vary in size depending on the number of relation tuples in $E_k$, we generate sets of varying sizes \footnote{In this work we use sizes 200, 500, 1K, 5K, 10K and 50K.} of ``positive'' random pairs $\{\langle n_i\ n_j\ 1\rangle \}$, where $n_i$ and $n_j$ are randomly sampled from $N_k$. We refer to these datasets as $\delta_{\mathrm{rand},x}^+$ where $x$ is the number of positive pairs generated.

 We need both positive and negative examples in the datasets to train a model. For ease of training, we aim to \textbf{generate balanced datasets} with the same number of positive and negative examples. Randomly \textbf{generating negative pairs} based on the seed vocabulary or selecting positive pairs of a different relation type is not optimal because models can learn to identify words/concepts associated to the relation rather than the relation itself~\cite{levy2015doEmbedsLearnRels}. Instead, we apply \emph{negative switching}, whereby positive pairs are switched based on the subject and object vocabularies for the relation, i.e. from positive pairs $\langle n_i\ n_j\ 1\rangle$ and $\langle n_k\ n_l\ 1\rangle$; thus subject vocabulary $\{n_i, n_k\}$ and object vocabulary $\{n_j, n_l\}$, it is possible to generate negative pairs $\langle n_k\ n_j\ 0\rangle$ and $\langle n_i\ n_l\ 0\rangle$. Albeit practical, this approach does not take into account transitive relations and assumes a closed world. Also, for relations with unbalanced subject-object vocabularies, it may be impossible to generate sufficient negative examples; in such cases we fall back to selecting pairs from other relation types or generating random pairs.

 \paragraph{Model training}
 We next use the generated datasets to train binary classification machine learning models for each studied embedding space $s \in S'_k$, we refer to the resulting trained models as $m^{\delta,s,t}$, where $\delta$ is the generated dataset and $t$ is the model type. We do not specify a type of machine-learning model to be used; however, since the input and output are very simple, we expect fully-connected neural networks to be suitable and use these in our experiments.

 Each sample $\langle n_i n_j l \rangle$ in the dataset is converted into an \emph{input vector for the model} by combining the embeddings for the subject and object arguments; i.e. $F_s(n_i) \odot F_s(n_j)$. In this work we use vector concatenation for $\odot$, but other operations are also possible. However, to prevent models from simply learning which embeddings are associated with a relation type, during training we apply \emph{input perturbation}: for each batch, we generate a random vector $v \in \mathbb{R}^d$ and add it to both the subject and object embedding. Thus the final input to the model is $(F_s(n_i) \oplus v) \odot (F_x{n_j} \oplus v)$. This ensures that the difference between the embeddings is the same as for the original pair, but no individual embedding is seen twice. Also, since different embedding spaces have different scales, we adapt the amount of perturbation to each embedding space by scaling the random perturbation to be within the standard deviation of the embedding space.

 Furthermore, to verify that the generated dataset does not encode information about the relation, we also train models using a baseline random embedding space $s_{\mathrm{rand}} \in S_k$ for the seed vocabulary. Since the embeddings in $s_{\mathrm{rand}}$ are random, they cannot encode any information about $n \in V_{\mathrm{seed}}$ and their relations $r \in R_k$.

 Each generated dataset is split into training, (validation) and testing subsets; the latter is used to evaluate the trained model by calculating metrics: precision, recall, accuracy and f1. Since the performance of a model is affected by random initialization of its parameters, we propose to perform multiple runs, resulting in $m_i^{\delta,s,t}$ for $1 \leq i \leq \#_{\mathrm{runs}}$. This allows us to calculate the average and standard deviation for each of the collected metrics. Below, we use $\mu^\upsilon_{\delta,s}$ to refer to the average for metric $\upsilon$ for models trained on dataset $\delta$ and embedding space $s$. If $s$ is omitted, the average is taken over all models trained on $\delta$. Likewise for the standard deviation $\sigma^\upsilon$.

 \paragraph{Prediction results analysis}
 Once we have trained all our models, we can do some statistic analysis on the collected metrics. Since the dataset generation and model training are mostly automated, our main goal here is to identify any non-predictable relations. Using our baselines, we can discard results based on two reasons: (i) the generated dataset is biased and (ii) the learned model's results are not significant.

 First, we define \textbf{baseline ranges for prediction metrics}. For this, we use the metrics gathered for the generated random datasets $\delta_{\mathrm{rand,x}}$ and define the range thresholds as $\tau^{\upsilon_\mathrm{min}}_{\mathrm{biased}} = \mu^\upsilon_{\delta_{\mathtt{rand},x}} - 2 \sigma^\upsilon_{\delta{\mathrm{rand},x}}$ and $\tau^{\upsilon_\mathrm{max}}_{\mathrm{biased}} = \mu^\upsilon_{\delta_{\mathtt{rand},x}} + 2 \sigma^\upsilon_{\delta{\mathrm{rand},x}}$. Any metrics within these ranges could be due to chance with 95\% probability.

 We consider a dataset to be \textbf{biased} if models can perform well on them regardless of whether the embeddings used to train the model encode any information. Intuitively, these are datasets which are imbalanced in some way allowing the model to exploit this imbalance during prediction, but that do not reflect the knowledge encoded in the embeddings. To detect these, we look at model results for models trained on random embeddings (i.e. on models $m^{\delta_r,s_{\mathrm{rand}},t}$). We say that $\delta_r$ is biased if $\mu^\mathrm{f1}_{\delta_r, s_{\mathrm{rand}}}$ is outside of the  $[ \tau^{\mathrm{f1}_\mathrm{min}}_{\mathrm{biased}}, \tau^{\mathrm{f1}_\mathrm{max}}_{\mathrm{biased}} ]$ range. The rationale is that even with random embeddings, such models were able to perform outside of the 95\% baseline ranges.

 We consider a trained model $m^{\delta_r,s,t}$ to be \textbf{significant} if its predictions are statistically better than predictions made by $m^{\delta_r,s_\mathrm{rand},t}$. Intuitively, this indicates that the embedding space $s$ contributes information about relation $r$ with high probability. Formally, we say that $m^{\delta_r,s,t}$ is significant if $\mu^{\mathrm{f1}}_{\delta_r,s} - \mu^{\mathrm{f1}}_{\delta_r, s_{\mathrm{rand}}} > 2\mathrm{max}(\sigma^{\mathrm{f1}}_{\delta_r,s}, \sigma^{\mathrm{f1}}_{\delta_r, s_{\mathrm{rand}}})$.

\section{Measuring Lexico-Semantic Knowledge in Embeddings}
\label{sec:lexicoSemanticStudy}

We applied the method described above to two lexico-semantic KGs (WordNet and Sensigrafo) and several embedding spaces derived from different corpora and embedding algorithms. In this section we describe how we applied the methodology and present a summary of the results obtained. In the next section we discuss analyze the results.\footnote{The code we used is available at \url{https://github.com/rdenaux/embrelassess}}

\subsection{Knowedge Graphs and Relations}
\label{sec:KGsAndRels}

We applied our methodology to two KGs: (i) WordNet, a well known lexico-semantic semantic network and (ii) Sensigrafo, a proprietary lexical knowledge graph developed by Expert System. Although similar in structure and scope --both KGs aim to provide a sense lexicon per language and use hypernymy as the main relation between senses--, Sensigrafo has been developed independently and is tightly coupled to Expert System's text analysis pipeline, enabling state of the art word sense disambiguation with claimed 90\% accuracy.

\begin{table}
  \centering
  \small
  \begin{tabular}{lll}
 & Sensigrafo & WordNet\\
\hline
version & 14.2 & 3.0\\
words/lemmas & 400K & 155K\\
syn(set/con) & 300K & 118K\\
relations & 55 & 27\\
derived rel datasets & 149 & 27\\
pair types & lem2(lem,syn,POS) & lem2lem\\
 & syn2(syn,POS) & \\
\hline
  \end{tabular}
  \caption{Lexical Knowledge Graphs used.}\label{tab:kgs}
\end{table}

Table~\ref{tab:kgs} provides an overview and comparison for the KGs. WordNet does not have an explicit identifier for concepts; instead, it defines sets of synonyms, i.e. lemmas with a shared sense. Sensigrafo calls such sets syncons (synonym-concepts), since they refer to specific concepts and assigns a unique identifier to each. There are differences in terms of size, granularity of concepts, structure of the network (e.g. choice of central concepts) and types of relations. Finally, because Sensigrafo is developed and maintained as part of a text analytics pipeline, some of its features are tailored and biased towards supporting  functionality and domains required by Expert System customers. By comparison, being a community effort WordNet may benefit from a wider set of stakeholders.

\begin{table*}
  \centering
  \small
  \begin{tabular}{llllrrrrr}
type & name & KG & example pair & lem pairs & syn  pairs & obj:subj & voc tot\\
\hline
Lexical & also see & w & mild-temperate & 5800 &  & 1.25 & 2339\\
 & derivation & w & revoke-revocation & 118888 &  & 1.0 & 18094\\
 & pertainym & w & regretfully-sorry & 6516 &  & 1.37 & 6079\\
     & participle of & w & operating-operate & 81 &  & 2.13 & 94\\
   \noalign{\smallskip}
Hypernymy & hypernym* & w & cinnamon-spice & 110650 &  & 0.45 & 27048\\
 & sup/subnomen & s & ditto & 56413 & 33696 & 1.98 & 19768\\
 & super/subverbum & s & kick-move & 37660 & 9426 & 1.81 & 6630\\
 & instance hypernym* & w & gemini-constellation & 2358 &  & 0.44 & 1526\\
   \noalign{\smallskip}
Categorical & category domain* & w & fly-air, tort-law & 9116 &  & 0.13 & 4166\\
 & sensiDomain & s & antitrust case-commercial law & 28610 &  & 0.02 & 17636\\
 & usage domain* & w & squeeze-slang & 846 &  & 0.09 & 395\\
 & region domain* & w & legionnaire-france & 1349 &  & 0.12 & 546\\
 & noun cat & s & flora-natural object &  & 186797 & 0.004 & 50257\\
 & verb cat & s & belong-v. of generic state &  & 27388 & 0.01 & 11498\\
 & tag & s & Christmas Eve-calendar day &  & 31775 & 0.01 & 16825\\
   \noalign{\smallskip}
Meronymy & member meronym* & w & archipelago-island & 1315 &  & 1.51 & 1147\\
 & substance meronym* & w & brine-sodium & 369 &  & 0.98 & 378\\
 & part meronym* & w & aeroplane-wing & 6403 &  & 1.39 & 4054\\
 & omni/parsnomen & s & construction-roofing & 807 & 3110 & 1.22 & 597\\
   \noalign{\smallskip}
Synonymy & synonym & w & encourage-promote & 74822 &  & 1.0 & 22304\\
 & synonym & s & ditto & 69974 &  & 1.0 & 19130\\
   \noalign{\smallskip}
Concept Simil & similar & w & big-immense & 20464 &  & 1.0 & 6790\\
 & attribute & w & short-length, good-quality & 1718 &  & 1.0 & 859\\
 & cause & w & secure-fasten & 719 &  & 0.93 & 445\\
 & syncon-cause & s & ring-sound, fright-fear & 584 &  & 0.87 & 458\\
 & entailment & w & look-see, peak-go up & 1519 &  & 0.89 & 956\\
 & syncon-implication & s & overtake-compete & 1358 & 291 & 0.83 & 870\\
 & verb group & w & shift-change, keep-prevent & 4944 &  & 1.0 & 1193\\
 & syncon-corpus & s & find-strike & 97707 & 39644 & 1.11 & 10502\\
 & syncon-unification & s & ritual-rite, enclose-envelop & 6599 & 2392 & 0.99 & 3366\\
 & antonym & w & release-detain & 9310 &  & 1.0 & 4651\\
 & (adj,n,v)antonym & s & ditto & 4656 & 1728 & 1.03 & 1823\\
   \noalign{\smallskip}
Positional & s-adjective-class & s & nightmarish-account & 131853 & 33037 & 0.69 & 13972\\
 & adverb-(n,v,adj,adv) & s & below-criteria & 15598 & 3367 & 1.12 & 2509\\
 & verb-object & s & counter-illness & 108046 & 23163 & 1.24 & 9087\\
 & verb-subject & s & less-tension, plunge-index & 46949 & 10131 & 1.01 & 7126\\
   \noalign{\smallskip}
Prepositional & 12 internoun & s & meeting-colleague &  &  & 0.83 & 10748\\
 & 10 verb prep noun & s & break down-tear &  &  & 0.56 & 1848\\
 & verb prep verb & s & set-rise & 287 &  & 1.19 & 196\\
 & 4 adj prep noun & s & eligible-admission &  &  & 1.62 & 299\\
   \noalign{\smallskip}
Geographic & geography & s & Brussels-Brussels Capit. Reg. &  & 1555 & 0.23 & 1584\\
   \noalign{\smallskip}
Part-of-Speech & Noun & s & entity-GNoun &  & 45961 & n/a & 76138\\
 & ProperNoun & s & Underground Railway-GPNoun &  & 4134 & n/a & 8273\\
 & Verb & s & cleanse-GVerb &  & 11343 & n/a & 22691\\
 & Adjective & s & record-GAdjective &  & 12585 & n/a & 25175\\
 & Adverb & s & in my opinion-GAdverb &  & 2110 & n/a & 4225\\
   \end{tabular}
  \caption{Overview of lexico-semantic relations studied.}
  \label{tab:relsA}
\end{table*}

WordNet and Sensigrafo provide 27 and 55 types of relations respectively which we used to generate datasets. Table~\ref{tab:relsA} shows an overview of the relations extracted; since discussing the 88 relations individually is unwieldy, we group the relations into types following --and expanding-- a categorization described in the WordNet manual\footnote{\url{https://wordnet.princeton.edu/documentation/wninput5wn}}. It suggests a top-level distinction between \textbf{Lexical} --those that hold between words-- and \textbf{Semantic} --those that hold between concepts-- relations and lists relations \emph{also see}, \emph{antonym}, \emph{derivation}, \emph{participle} and \emph{pertainym} as lexical. However, \emph{antonym} is clearly based on the word's meaning. \emph{Derivation} and \emph{participle} seem to be truly lexical, while the others can also contain pairs which are semantically related, as shown by the example pairs in the table. Sensigrafo does not contain purely lexical relations.

The main relation type used by both KGs is \textbf{hypernymy}, which relates narrower to broader concepts (or instances to concepts). Although hypernymy is transitive, we only consider direct relations explicitly stated in the KG during the dataset generation.

\textbf{Categorical relations} are those that associate a concept with some category. When categories are concepts, they can be seen as an indirect hypernymy. WordNet has three types of categories called domains: \emph{category} --concept to a topic domain--, \emph{usage} --concepts to a type of use-- and \emph{region} --concepts to places. Sensigrafo has a set of core noun and verb concepts --called noun or verb categories, as well as tags--  which form the backbone of the hypernymy hierarchy. These category syncons tend to be quite abstract concepts, hence we have only considered these relations at the syn2syn level, as they do not have lemmas in the seed vocabulary. Sensigrafo also defines a list of about 400 domains (which are not syncons), which is similar to WordNet's \emph{category domain}.

\textbf{Meronymy} relates concepts in whole-member relations. WordNet distinguishes between \emph{membership} --part of group--, \emph{substance} --when something is made of substances that are orders of magnitude smaller than the whole-- and \emph{part} -- remaining cases. Sensigrafo does not distinguish between these cases and uses this type of relation sparingly.

\textbf{Synonymy} can be defined between lemmas in the same synset/-con and is a special case of \textbf{conceptual similarity}. The \emph{similarity} relation in WordNet captures pairs of similar adjectives. The \emph{attribute} relation in WordNet relates adjectives describing a value for a noun (similar, but more limited than the \emph{adjective-class} relation in Sensigrafo). The \emph{cause} relation in WordNet describes causality between verbs and is similar in Sensigrafo. The \emph{entailment} relation in WordNet describes an entailment between verbs, similarly to the \emph{syncon-implication} in Sensigrafo, which can also be applied to non-verbs. \emph{Verb group} in WordNet groups similar verbs. \emph{Syncon-unification} is assigned by linguists to syncon pairs that could be merged as a single syncon. \emph{Antonym} describes concepts with opposite meanings.

\textbf{Positional} relations encode co-locations of concepts and are only provided by Sensigrafo. Two main subtypes: one relates adjectives or adverbs to other concepts while the second relates verbs to concepts that appear as the verb subject or object.

Sensigrafo provides about two dozens of \textbf{prepositional} relations between concepts. Such relations are of the type POS+preposition-POS. E.g. the pair \emph{rival-titleholder} has relation \emph{noun+to-noun}. Sensigrafo also encodes \textbf{geographic} relations between places. However, since the seed vocabulary did not contain many place names, we could not generate a dataset between lemmas.

Finally, we also generated datasets relating syncons to their \textbf{part-of-speech} (POS) besides the various $\delta_{\mathrm{random},x}$ as explained in Section~\ref{sec:genericMethod} (not included in Table~\ref{tab:relsA} ).

We generated a total of 176 datasets, based on the 88 relations and a seed vocabulary consisting of 76K concepts and 71K lemmas (roughly corresponding to our smallest embedding space, which was trained on a disambiguated version of the English United Nations corpus~\cite{ziemski2016united}). For WordNet, we only generated datasets between lemmas. For Sensigrafo we also generated datasets between words, word/concepts and concepts.

\subsection{Embeddings and Corpora}

The word and concept embeddings studied were derived from 6 algorithms. Three provided embeddings based on sequences of either words or syncons: GloVe~\cite{pennington2014glove}, FastText~\cite{Bojanowski2017}, Swivel~\cite{Shazeer2016}. Two provided joint word and concept embeddings: Vecsigrafo~\cite{denaux2017vecsigrafo} based on a disambiguated corpus and HolE~\cite{nickel2016HolE}, which directly generates embeddings from a KG (not from a corpus) and thus serves as a reference point. The embeddings were either publicly available or provided to us by \cite{denaux2017vecsigrafo}. We also generated random embeddings $s_{\mathrm{rand}}$ for the seed vocabulary.

The embeddings were trained on three different corpora, which we chose to study whether relation prediction capacity varies depending on the corpus size: the English United Nations corpus\cite{ziemski2016united} (517M tokens), the English Wikipedia (just under 3B tokens) and Common Crawl (around 840B tokens). Although we aimed at using embeddings with 300 dimensions, in a few cases we had to diverge as the embeddings were only available with other dimensions.

\subsection{Training and Results}
We generated models based on fully connected neural networks (NN) with either 2 or 3 hidden layers. Initially we also generated models with logistic regression, which consistently underperformed. The 2 and 3 layer NNs generally produced similar results, suggesting they converge for the given datasets. For the standard case of embedding dimension 300, the input to the net is a vector of 600 dimensions, followed by hidden layers of 750 and 400 nodes for the NN2; and (750, 500, 250) for the NN3 (similar architectures were defined for single embedding relations such as syn2POS). The output layer has 2 nodes and uses 1-hot-encoding to encode positive or negative examples. For regularization, we use dropout between layers with value 0.5. We used a heuristic rule that varies the number of epochs to train a dataset depending on its size. Datasets with $<300$ positive examples are trained for 48 epochs, those with $<5K$ for 24, $<30K$ for 12 and large datasets are trained for 6 epochs. We use the Adam optimizer with learning rate $1^{-5}$ and the cross entropy loss function. A scheduler reduces the learning rate on plateau. All of these hyper-parameters where derived through trial and error with a few sample relation datasets and kept constant to automatically train models without manual inspection. We used a random 90, 5, 5 split for training, validation and test from the input dataset. For WordNet relations we used a mixture of NN2 and NN3 models and trained each model either 3 or 5 times. As training the models is the main bottleneck of our approach, for Sensigrafo relations we only trained NN3 models, training each model 3 times.

We trained a total of 10,560 models, resulting in 1,596 evaluation metrics averaged over $n$ runs: 126 for the random relation datasets, 149 trained on random embeddings, 1,065 for Sensigrafo relations and 268 for WordNet relations. Each run resulted in metrics for relation prediction on unseen pairs during training.

\section{Analysis and Discussion of Results}
\label{sec:analysis}

\begin{table*}
  \centering
  \small
  \begin{tabular}{ll @{\hskip 0.25in}rr @{\hspace*{5mm}}rrr@{\hskip 0.2in}rr @{\hskip 0.2in}rr}
    \hline
    &  & \multicolumn{2}{c}{}   & \multicolumn{5}{c}{models} & \multicolumn{2}{c}{metrics}\\
    \cline{5-9}
    dataset rel & KG & \multicolumn{2}{c}{datasets}   & & \multicolumn{2}{c}{non-predictable} & \multicolumn{2}{c}{``predictable''} & absolute & relative \\
    \cline{3-4}\cline{6-7
 &  & \# & \# biased & \# & \% biased & \% not signif. & \textbf{\% better} & \% worse & $\mu_{\mathrm{f1}}$ & $\Delta \mu_{\mathrm{f1}}$\\
\hline
all & both & 156 & 38 & 1281 & \textbf{28.6} & \textbf{46.5} & \textbf{23.9} & \textbf{0.9} & 0.687 & 0.173 \\[3pt]
all & wn & 19 & 6 & 216 & 27.8 & 39.4 & 30.1 & 2.8 & 0.684 & 0.149\\
all & sensi & 137 & 32 & 1065 & 28.8 & 48. & 22.6 & 0.6 & 0.688 & 0.181\\[4pt]
concept & sensi & 44 & 10 & 293 & 22.9 & \textbf{73.7} & \textbf{3.4} & 0. & \textbf{0.870} & \textbf{0.608}\\
word/concept & sensi & 43 & 3 & 172 & \textbf{7.} & 43.6 & \textbf{48.8} & 0.6 & 0.664 & 0.162\\
word & sensi & 50 & 19 & 600 & \textbf{38.} & 36.7 & 24.5 & 0.8 & 0.690 & 0.162\\[4pt]
lexical\(_{\text{w}}\) & wn & 4 & 0 & 32 & 0. & \textbf{75.} & 25. & 0. & 0.652 & 0.167\\[2pt]
hypernym\(_{\text{c}}\) & sensi & 2 & 1 & 14 & 50. & 42.9 & 7.1 & 0. & \textbf{0.916} & \textbf{0.696}\\
hypernym\(_{\text{w/c}}\) & sensi & 2 & 1 & 8 & 50. & 0. & \emph{50.} & 0. & 0.699 & 0.148\\
hypernym\(_{\text{w}}\) & sensi & 2 & 1 & 24 & 50. & 0. & \emph{50.} & 0. & 0.713 & 0.143\\
hypernym\(_{\text{w}}\) & wn & 2 & 1 & 24 & 50. & 4.2 & 37.5 & \textbf{8.3} & \textbf{0.756} & 0.106\\[2pt]

categ\(_{\text{c}}\) & sensi & 3 & 1 & 21 & 33.3 & 52.4 & 14.3 & 0. & \textbf{0.891} & \textbf{0.671}\\
categ\(_{\text{w/c}}\) & sensi & 4 & 1 & 16 & 25. & 12.5 & \textbf{62.5} & 0. & 0.744 & 0.263\\
categ\(_{\text{w}}\) & sensi & 1 & 1 & 12 & \textbf{100.} & 0. & 0. & 0. &  & \\
categ\(_{\text{w}}\) & wn & 4 & 4 & 40 & \textbf{100.} & 0. & 0. & 0. &  & \\[2pt]

meronym\(_{\text{c}}\) & sensi & 2 & 0 & 14 & 0. & \textbf{92.9} & 7.1 & 0. & 0.667 & 0.020\\
meronym\(_{\text{w/c}}\) & sensi & 1 & 0 & 4 & 0. & \textbf{50.} & \emph{50.} & 0. & 0.664 & 0.296\\
meronym\(_{\text{w}}\) & sensi & 1 & 0 & 12 & 0. & \textbf{91.7} & 0. & \textbf{8.3} &  & \\
meronym\(_{\text{w}}\) & wn & 3 & 0 & 48 & 0. & \textbf{66.7} & 31.2 & 2.1 & 0.708 & 0.166\\[2pt]

synon\(_{\text{c}}\) & sensi & 0 & 0 & 0 & 0. & 0. & 0. & 0.&  & \\
synon\(_{\text{w/c}}\) & sensi & 1 & 0 & 4 & 0. & 25. & \textbf{75.} & 0. & \textbf{0.804} & 0.249\\
synon\(_{\text{w}}\) & sensi & 1 & 0 & 12 & 0. & 8.3 & \textbf{91.7} & 0.  & 0.677 & 0.114\\
synon\(_{\text{w}}\) & wn & 1 & 0 & 8 & 0. & 0. & \emph{100.} & 0. & 0.680 & 0.135\\[2pt]

simil\(_{\text{c}}\) & sensi & 6 & 1 & 42 & 16.7 & 81. & 2.4 & 0. & \textbf{0.909} & \textbf{0.688}\\
simil\(_{\text{w/c}}\) & sensi & 7 & 0 & 28 & 0. & 57.1 & 42.9 & 0.  & 0.624 & 0.210 \\
simil\(_{\text{w}}\) & sensi & 7 & 1 & 84 & 14.3 & 38.1 & \textbf{47.6} & 0. & 0.628 & 0.128\\
simil\(_{\text{w}}\) & wn & 5 & 1 & 64 & 12.5 & 43.8 & 39.1 & \textbf{4.7} & 0.655 & 0.153\\[2pt]

position\(_{\text{c}}\) & sensi & 6 & 1 & 42 & 16.7 & \textbf{78.6} & 4.8 & 0. & \textbf{0.886} & 0.667\\
position\(_{\text{w/c}}\) & sensi & 6 & 1 & 24 & 16.7 & 37.5 & \textbf{45.8} & 0. & 0.721 & 0.120\\
position\(_{\text{w}}\) & sensi & 7 & 5 & 84 & \textbf{71.4} & 21.4 & 6. & 1.2 & 0.703 & 0.069\\[2pt]

prepos\(_{\text{c}}\) & sensi & 19 & 4 & 133 & 21.1 & \textbf{78.9} & 0. & 0. &  & \\
prepos\(_{\text{w/c}}\) & sensi & 22 & 0 & 88 & 0. & 51.1 & \textbf{47.7} & 1.1  & 0.629 & 0.124\\
prepos\(_{\text{w}}\) & sensi & 27 & 11 & 324 & 40.7 & 48.5 & 9.9 & 0.9 & 0.677 & 0.097\\[2pt]

POS\(_{\text{c}}\) & sensi & 5 & 1 & 20 & 20. & \textbf{70.} & 10. & 0. & \textbf{0.882} & \textbf{0.665}\\
POS\(_{\text{w/c}}\) & sensi & 0 & 0 & 0 & 0. & 0. & 0. & 0. &  & \\
POS\(_{\text{w}}\) & sensi & 4 & 0 & 48 & 0. & 2.1 & \textbf{97.9} & 0. & \textbf{0.746} & 0.262\\
\hline
\end{tabular}
  \caption{Overview of results.}
  \label{tab:resultsOverview}
\end{table*}

\paragraph{Biased datasets and non-significant models}
We apply the prediction result analysis  described in Section~\ref{sec:genericMethod}: based on 12 random datasets $\delta_{\mathrm{rand},x}$ and 126 prediction results we obtained $\mu^{\mathrm{f1}}_{\delta_{\mathrm{rand}}}=0.41$ and  $\sigma^{\mathrm{f1}}_{\delta_{\mathrm{rand}}}=0.12$, resulting in a baseline range for prediction metrics (i.e. $[ \tau^{\mathrm{f1}_\mathrm{min}}_{\mathrm{biased}}, \tau^{\mathrm{f1}_\mathrm{max}}_{\mathrm{biased}} ]$) of $[0.16, 0.65]$.

Using this range, we identify 38 of 156 datasets as being \emph{biased}. With other words, even though we took care not to generate non-biased datasets by using negative switching, a little more than a quarter of the datasets generated in this way contains clues about the relation. Interestingly, word-pair datasets are much more likely to be biased (38\%) while word/concept pairs are unlikely to so (only 7\%) (see Table~\ref{tab:resultsOverview}). These result suggest that using KGs as a \emph{silver standard}~\cite{Paulheim2017KGRefinement} is of limited use for word-pair prediction, but is suitable when linking words to concepts.

We see that using a KG to build datasets for certain kind of relations is very difficult. This is the case for word pair datasets for categorical relations (all of the 52 generated datasets were biased); this is likely due to such relations being highly unbalanced with only a few words in the object position, further compounded by ambiguity of words compared to concepts. We think word ambiguity also plays a role in the difficulty producing non-biased datasets for positional relations.

On the bright side, about 71\% of the models were trained on a non-biased dataset and we use these to identify (non)significant models. Overall, we found that 46.5\% of the trained models using pre-trained embeddings were not statistically significant different from a baseline using random embeddings. About 1\% of the models performed worse than the baseline. These were typically small datasets and relations that were hard to learn. In any case this 1\% is below the 2.5\% that can be expected by using the $2\sigma$ threshold for significance.

Conversely, only about 24\% of the models trained using pre-trained embeddings significantly outperformed the baseline. For word/concept pair datasets this percentage jumps to almost 49\% while for concept pairs it plummets to only 3.4\%. Since most of the pre-trained embeddings we are using are corpus-based, this shows that such embeddings have trouble capturing the relations at the purely conceptual level. At the same time, these models are relatively good at relating words to concepts using semantic relations. In the sections below we only discuss results for the models that outperformed their baseline.

\subsection{Relation Types in Embeddings}

In the last two columns of Table~\ref{tab:resultsOverview} we see the absolute and relative f1 measures for the models that significantly outperformed the random baseline. We see that most of the results with an average f1 score higher than 0.8 are for relations between concepts; these were all achieved by training models on the HolE embeddings. Swivel embeddings also obtained good results for predicting categorical relation (0.846 f1, but this result just cleared the $2\sigma$ significance threshold). One model trained on FastText embeddings achieved (0.666 f1 on a meronymy relation). This confirms that \emph{the studied corpus-based embeddings are not capable of capturing semantic relations at the concept level}.

For the models trained on word/concept or word pairs, we can go through the relation types. Absolute prediction accuracy for \textbf{lexical relations} is poor, the Vecsigrafo and GloVe are the only embeddings that can predict pertainym relations with over 0.7 f1 score. Prediction of \textbf{hypernym} relations is good; the best performing model is trained on HolE on a word/concept dataset (0.797 f1); from the corpus-based embeddings, FastText performs best (0.766). Prediction of \textbf{categorical} relations (only for word/concept pairs) is quite good; corpus-based embeddings using Vecsigrafo perform similarly to HolE (around 0.82 f1) for nouns; corpus-based embeddings have trouble with categorical relations between verbs, with performance dropping below 0.6 compared to 0.8 for HolE. For \textbf{meronymy} relations, the datasets generated on Sensigrafo were rather small, resulting in f1 performances around 0.69 with Vecsigrafo embeddings; WordNet-derived datasets produce good performance with FastText, GloVe and Vecsigrafo having scores between 0.75 and 0.81 for the part-meronym relation; substance- and member-meronym relations accuracy drops below 0.7. Performance for \textbf{synonym} relations depends strongly on the used embedding; the best performers are HolE (0.93 f1), Vecsigrafo (0.79) and FastText (0.75) while Swivel performs poorly. Prediction of \textbf{similarity} relations is mediocre; the best embeddings are HolE (0.81), followed by GloVe and FastText (around 0.74); depending on the relation, prediction f1 can drop to 0.6. For \textbf{positional} relations, Vecsigrafo (0.76) outperforms even HolE (0.74) with FastText and Glove trailing (0.7); this is similar for \textbf{prepositional} relations. For \textbf{part-of-speech}, FastText and GloVe perform well (f1 above 0.8).

\subsection{Impact of embedding type and corpus size}

Unsurprisingly, the larger the corpus the better the overall results for Sensigrafo datasets. Average f1 scores for embeddings trained on the UN corpus, Wikipedia and Common Crawl were 0.66, 0.69 and 0.73. For WordNet datasets (only lem2lem), the scores are similar: 0.68, 0.72 and 0.71. Thus, although increasing the corpus size helps, for many types of relations, the gain from training on a very large corpus is relatively small.

Table~\ref{tab:algoByPairTypes} summarizes the performance of the different embedding learning algorithms for predicting relations grouped by the pair type. For word pair prediction, FastText outperforms other embedding learning algorithms, including HolE (although this difference is not major). HolE excels at predicting relations between senses, but its performance decreases as lemmas are introduced, since it cannot disambiguate between the senses. Vecsigrafo and GloVe both are not far behind the performance of FastText and HolE, but produce significant predictions for more relations than FastText and HolE. Standard Swivel with words lags behind, especially in the number of relations that it can predict.

\begin{table}
  \centering
  \small
\begin{tabular}{rlrrl}
datasets & algo & F1\(_{\text{avg}}\) & F1\(_{\text{std}}\) & pair type\\
\hline
\textbf{8} & HolE & \textbf{0.90} & 0.02 & concept\\
1 & Swivel & 0.85 & 0.0 & concept\\
1 & FastText & 0.67 & 0.0 & concept\\
  \noalign{\smallskip}
26 & HolE & \textbf{0.67} & 0.09 & word/concept\\
48 & Vecsigrafo & 0.64 & 0.08 & word/concept\\
  \noalign{\smallskip}
38 & FastText & \textbf{0.74} & 0.08 & word\\
12 & HolE & 0.72 & 0.09 & word\\
43 & Vecsigrafo & 0.68 & 0.07 & word\\
43 & GloVe & 0.68 & 0.08 & word\\
21 & Swivel & 0.66 & 0.06 & word\\
\end{tabular}
  \caption{Average F1 scores for predicting relations with different pair types and embeddings.}\label{tab:algoByPairTypes}
\end{table}

\subsubsection{Impact of joint word-concept learning}
Vecsigrafo co-trains word and concept embeddings. Table~\ref{tab:algoByPairTypes} shows that compared to Swivel, this co-training improves, for word pairs, both the number of relations that can be predicted (double the number) and the average f1 score (we assume that the additional predicted relations push the average score down). Similarly, we see that compared to word pairs, Vecsigrafo is able to produce predictions for more relations when considering word/concept pairs. As discussed above, corpus-based, joint word-concept training does not seem to capture relations between at the concept level, suggesting there is room to improve such algorithms.

\section{Conclusion}
This paper presented a methodology for studying whether embeddings capture relations as well as KGs and applied it to study lexico-semantic relations between words and concepts. The results show that for a few relations, word embeddings can outperform embeddings derived directly from KGs. Also, corpus-based embeddings fail to capture relations at the concept level. However, for most relation types, embeddings only can predict relations with an accuracy under 0.7. Our results provide evidence that correct capture of relations should happen at the concept level and may not be achievable with high accuracy at the word level. As future work, we want to apply our method to study contextual embeddings.

\section*{Acknowledgments}
The research reported in this paper is supported by the EU Horizon 2020 programme, under grants European Language Grid-825627 and Co-inform-770302.

\bibliographystyle{named}


\section*{}
\vspace{-1cm}

\footnotetext{\textit{$^{a}$~Ural Federal University
620002, 19 Mira street, Ekaterinburg, Russia. Tel: +7 343 375 4444; E-mail: contact@urfu.ru}}
\footnotetext{\textit{$^{b}$~College of Science, Institute of Materials Physics and Chemistry, Nanjing Forestry University, Nanjing 210037, PR China }}

\section{INTRODUCTION}
Monolayer VSe$_2$ is the one of the most intriguing members of the family of two-dimensional (2D) transition-metal dichalcogenides. This material attracts a special interest of the scientific community due to several recent discoveries, including in-plane piezoelectricity \cite{1}, a pseudogap with Fermi arc \cite{2} at temperatures above the charge density wave transition (~220 K for the monolayer \cite{3}), and especially the existence of ferromagnetism in 2D system\cite{4,5,6,7,8,9,10,11}. Experimental results are rather contradictory. A strong room-temperature ferromagnetism with a huge magnetic moment per formula unit has been reported for monolayer VSe$_2$ epitaxially grown on graphite \cite{4}. A local magnetic phase contrast has also been observed by magnetic force microscopy at the room temperature at the edges of VSe$_2$ flakes exfoliated from a three-dimensional crystal. \cite{14} XMCD measurements evidence a spin-frustrated magnetic structure in VSe$_2$ on graphite. \cite{xmcd} Paramagnetism of bulk VSe$_2$ \cite{para1,para2} makes these observations more intriguing. Another situation was reported for the monolayers grown on bilayer graphene/silicon carbide substrate. In both works the absence of exchange splitting of the vanadium $3d$ bands observed in angle-resolved photoemission spectroscopy experiments was reported. This result contradicts to other studies that revealed a magnetization value not higher than ~5 $\mu_B$.\cite{12,13} Based on these results we can conclude that the influence of the substrate is important for description of the magnetic properties of these materials. Theoretical models have been developed to account for the above discrepant observations \cite{4,12,14,16}. These works mainly focused on the band structure and magnetic moments on vanadium sites. It has been proposed that the presence of charge density waves could cause the quenching of monolayer ferromagnetism due to the band gap opening induced by Peierls distortion \cite{15}. Phonon spectra of several VSe$_2$ and similar systems
were also considered theoretically \cite{22,23}. This modeling motivates us to study interplay between magnetism and structural phase transitions in VSe$_2$. Additionally, there is a plethora of works demonstrating a relationship between the symmetry, electronic structure and magnetic properties in transitional metal compounds \cite{BaCu, DMI1, DMI2,24}.

The VSe$_2$ crystal is formed from separate layers along the c-axis direction. Two main phases for this material were predicted to be stable: the H phase characterized by Se stacked over each other and the T phase with one layers of Se rotated by $60^\circ$ around axis normal to the plane of layer. \cite{16} Atomic structures of the VSe$_2$ monolayer in both H and T phases are shown in Fig.1. Surprisingly, the reported binding energies for different configurations are almost the same despite the colossal difference in magnetic properties and electronic structure (Fig.\ref{fig1}). \cite{16} This finding additionally motivates us to examine various aspects of structural phase transitions in bulk, few-layer and monolayer of VSe$_2$.
\newpage
\begin{figure}[t]
    \centering
    \includegraphics[width=1\columnwidth]{Fig1.png}
    \caption{Atomic structure of 2D VSe$_2$ monolayer (top and side view) in H phase (a) and in T phase (b). Vanadium atoms are denoted with red circles, upper and bottom selenium layers denoted with light green and dark green circles, respectively. (c) and (d) panels represent the corresponding spin-polarized band structures. Red lines correspond to spin up states and black ones to spin down, the Fermi level corresponds to 0 eV.}
    \label{fig1}
\end{figure}
\section{Computational method and model}
Electronic properties of the VSe$_2$ system were
simulated within Density Functional Theory (DFT) framework using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional \cite{17} as implemented in the Vienna ab-initio simulation
package (VASP) \cite{19,20} with a plane-wave basis set. This approach gave reliable results for other systems similar to VSe$_2$ \cite{25}. Also we include van der Waals interaction using the method of Grimme (DFT (PBE)-D2) \cite{vdW}. Taking into account London
dispersion forces is essential for few-layer VSe$_2$ (see Table
\ref{tab1} and discussion in section 3.5).

The calculation parameters were chosen as follows. The energy cutoff
equals to 400 eV and the energy convergence criteria is
$10^{-6}$ eV. For the Brillouin zone integration a $10\times10\times1$
gamma centered grid was used for layered structures and $8\times8\times8$ for bulk structures. A vacuum space more than 10 \AA \, in the vertical $z$ direction
was introduced for layered structures. The technical parameters are similar to those used in the recent studies of phase stability in layered systems.\cite{Ersan, Kaltsas}

The optimized atomic positions for T-phase and lattice parameters $a=b=3.31$ \AA \, and $c=6.20$ \AA \, are in good agreement with experiment \cite{21}.
In particular, the corresponding interlayer distance in bulk VSe$_2$ is $3.04$ \AA. The calculated band structures of VSe$_2$  monolayer in the T and H phases are in good agreement with previous works. \cite{16} The calculated magnetic moment of 0.68 $\mu_B$ for initial configuration without rotation of the selenium atoms also agrees with results of the previous work \cite{18}.

To investigate the transition between H and T phases we performed self-consistent calculations of electronic structure and total energies in transitional points between these phases. For this purpose, we rotate either one Se atom or all selenium atoms belonging to the upper layer of VSe$_2$ in supercell as schematically shown in Fig.\ref{fig2}. To trace the changes in electronic structure and magnetic properties the calculations for configurations with a 10$^{\circ}$ rotation step were performed. Generally, the rotation can be realized within two models. The first one is to move Se in plane from initial to final point (Fig.\ref{fig2} a and c). The second one is to fix the constant V-Se distance for all intermediate steps, which produces an elevation of selenium atoms above the plane at intermediate steps of the migration (Fig.\ref{fig2} b and d). We will refer these rotation models as in-plane and arc rotation schemes, respectively. All the calculations were performed for the ferromagnetic ordering of the spins of vanadium atoms.
\begin{figure}[h]
\includegraphics[width=0.875\columnwidth]{Fig2.png}
\caption{Schematic visualization of the plane (a,c) and arc (b,d) types of the Se atoms rotation. (a,b) and (c,d) panels correspond to side and top views, respectively. Initial and final positions of Se are presented with orange and green circles, respectively. Intermediate configurations of selenium atoms obtained with the 20$^{\circ}$ step are denoted with light blue circles.}
    \label{fig2}
\end{figure}

\section{Results and discussion}
\subsection{Rotation of single Se atom}
At the first step of our study we have simulated the motion of the single Se atom in the monolayer (see Fig. \ref{fig3}). For simplicity, we considered an in-plane migration of the atom. Results of the calculations (Fig. 3) evidence a gradual increasing of the total energy of the system during the all processes of the rotation with maximal value at final point. The cause of the large magnitude of the energies and instability of the final configuration is in decreasing of the distance between moved and rigid Se atoms to the value of 1.92 \AA. Thus we can conclude that the model of the single Se atom rotation is unrealistic and transition between T and H phases may be realized only with distortion of the whole selenium layer. Further we will consider only this kind of the structural phase transitions. The values of the magnetic moments calculated for intermediate configurations (Fig. \ref{fig3}) support our initial guess that the structural transition between the phases affects magnetic properties of VSe$_2$. Note that a deviation of the selenium atoms from equilibrium positions on small angles (less than 10$^{\circ}$) requires much smaller energies of about 0.32 eV and, therefore, should be taken into account for a realistic description of the atomic structure of VSe$_2$ at the room temperature.

\begin{figure}[ht]
    \centering
    \includegraphics[width=1\columnwidth]{Fig3.png}
    \caption{Evolution of the total energy (a) and magnetic moment (b) during in-plane rotation of single Se atom. (c) and (d) panels visualize the initial and final atomic structures. Light and dark green circles denote upper and bottom selenium layers, respectively.}
    \label{fig3}
\end{figure}
\subsection{Rotation of the whole Se sheet in the VSe$_2$ monolayer}
 Having considered the results concerning the migration of the single Se atom we are in a position to analyze the case of whole upper Se-layer rotation, which will provide a better understanding of the transition between H and T phases.
 \begin{figure}[ht]
    \centering
    \includegraphics[width=1\columnwidth]{Fig4.png}
    \caption{Evolution of the total energy (a) and magnetic moment (b) during rotation of whole upper Se layer of VSe$_2$ monolayer within arc model.}
    \label{fig4}
\end{figure}

 The performed simulations for 3$\times$3 supercell with constant V-Se distances when Se atoms elevate from initial and starting positions (Fig. \ref{fig2}) revealed the energy barrier of 0.60 eV that is smaller than that observed in the case of the in-plane rotation (Fig.\ref{figS1}a in SI). Thus further we will consider only this type of the Se atom migration. To evaluate the temperature required to overcome this barrier one should establish a relation between the calculated energies of the process and temperature of the reactions. We have addressed this question in our previous work \cite{go} and found that the barrier values of about 0.50 eV and 1.20 eV correspond to the room temperature and 200 $^{\circ}$C, respectively. Thus, the energy barrier of 0.60 eV can be overcome already at the temperatures about 40 $^\circ$C.

 Four conclusions could be drawn from these results. (i) There is a possibility of the structural phase transition in previously studied VSe$_2$ samples during measurements. (ii) For development of devices based on VSe$_2$ and similar monolayer systems one should take into account possibility of the structural phase transitions caused by the heating of the devices during work. Such a transition can significantly affect the work of the device due to difference in electronic structures of different  phases (see Fig. \ref{fig1},also changes in band structure Fig.\ref{figS3} in SI). (iii) One can use VSe$_2$ and similar systems as temperature detectors. (iv) According to our results there is a low-energy cost to deviate the selenium atoms belonging to one layer on a small angle from the equilibrium positions. It means that one needs to account this for a realistic interpretation of the experimental data.

Moderate temperature of the transition between different structural phases requires an examination of the electronic structure and magnetic properties at intermediate steps of the structural phase transition. The obtained calculations results demonstrate that in the case of the ferromagnetic ground state the values of magnetic moments change gradually with small step of 10$^\circ$ of the rotation of Se layer. From Fig.\ref{fig4} one can see that at 30$^\circ$ the magnetic moment has the maximal value of 1.05$\mu_{B}$, which is about two times larger than that in the initial configuration. According to the calculated occupation matrices such a magnetic moment change is mainly related to the contributions of $xy$ and $x^2-y^2$ orbitals of vanadium atoms (see Fig. \ref{fig5}). Since the total occupation (spin-up + spin-down) of the different orbitals remain almost the same, the orbital magnetic moment values change is fully connected with a redistribution of the electrons between different spin channels due to change of the hybridization between V and Se.
\begin{figure}[ht]
    \centering
    \includegraphics[width=1\columnwidth]{Fig5.png}
    \caption{Partial densities of states calculated for VSe$_2$ monolayer in the ferromagnetic configuration. The arc rotation scheme  with the $20^{\circ}$ step was used. Left and right panels correspond to ($d_{xy}$, $d_{yz}$, $d_{xz}$) and ($d_{x^2-y^2}$, $d_{z^2}$) sets of states, respectively.   }
    \label{fig5}
\end{figure}

In the case of an antiferromagnetic configuration the situation is more complicated. First of all, the magnetic lattice of VSe$_2$ is frustrated one, which is in agreement with experimental observations.\cite{xmcd}  This means that within a mean-field DFT approach we cannot define an antiferromagnetic collinear-type order corresponding to minimum of the magnetic interaction energy for all V-V bonds, simultaneously. The second complication follows from the fact that the system in question is metal. It means that the magnetization of individual vanadium atom can be very sensitive to the orientation of the neighbouring magnetic moments \cite{FeGe}. Indeed, our DFT simulations of the VSe$_2$ supercell with antiferromagnetic ordering have revealed a strong suppression of the magnetic moment values of some vanadium atoms in the supercell. In addition, we observe that the details of the magnetic moments suppression strongly depend on the size of the supercell. In this complex situation some information on magnetic couplings in the VSe$_2$ system could be extracted by using the theory of infinitesimal spin rotations approximation \cite{FeGe, Liechtenstein}. However, the magnetic couplings calculated in this way can be used for analysis only in the vicinity of the ferromagnetic configuration.

The values of the magnetic moments in AFM phase can be stabilized by inclusion of the on-site Coulomb interaction as can be done with DFT+$U$ approach. However, the using of the DFT+$U$ approach in the case of VSe$_2$ is questionable, since the experimental ARPES spectra are in good agreement with GGA band structure as it was shown in Refs.\cite{xmcd,arp1,12}. At the same time the inclusion of the Hubbard $U$ leads to considerable changes in the band structure.

Thus, the energy difference between AFM and FM solutions for VSe$_2$ simulated with GGA does not allow us to construct a comprehensive magnetic model and estimate the corresponding magnetic interactions between vanadium atoms. Nevertheless, the results of these calculations evidence that despite the changes of electronic structure at intermediate steps the ferromagnetic configuration remains significantly energetically favorable in all the cases (Fig.\ref{fig6}). Thus the possible structural distortions in VSe$_2$ will not provide a suppression of ferromagnetism. Our calculations demonstrate that possible transition from experimentally observed T phase toward H phase should provide an enhancement of ferromagnetic interactions and increasing of magnetic moment. To simulate the experimentally observed paramagnetic state of bulk VSe$_2$ \cite{para1,para2} one can use a dynamical mean-field theory.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.92\columnwidth]{Fig6.png}
    \caption{Difference of AFM and FM state total energies calculated in arc and plane rotation schemes for 3x3 supercell of  VSe$_2$ monolayer.}
    \label{fig6}
\end{figure}

\subsection{Structural phase transition in bulk VSe$_2$}
There are two main differences in the energetics of the structural phases transitions in bulk and monolayer VSe$_2$. The first one is the almost the same value for the energies of the motion of Se layer within both rotation models (Fig.\ref{fig7} and Fig.\ref{figS1}c in SI). The second one is increasing of the migration barrier (see Fig. \ref{fig7}). Both are related to the van der Waals interactions between the layers in the bulk VSe$_2$. The analysis of the calculated partial density of states in this case leads to similar conclusions as above (see Fig.\ref{figS2} in SI)
\begin{figure}[ht]
    \includegraphics[width=1\columnwidth]{Fig7.png}
    \caption{Total energy (a) and magnetic moment (b) as functions of the rotation angles. The simulation results were obtained for bulk VSe$_2$ within the arc rotation model.}
    \label{fig7}
\end{figure}

In the case of the rotation of the Se atoms belonging to one layer with constant V-Se distance at intermediate steps of migration, the initial distance of 3.63 \AA \, between rotated and fixed selenium layers decreases by 0.54 \AA. This deviation from the optimal interlayer distance provides an increasing of the energy barrier (see also changes in band structure Fig.\ref{figS4} in SI).  The value of the energy barrier is corresponding with stability of the structural ground state in bulk crystal up to the temperatures above 100$^{\circ}$C. Note that in contrast to monolayer case the structural ground state of bulk VSe$_2$ is T configuration with  ferromagnetic orientation of magnetic moments.

\subsection{Structural phase transition in bi- and trilayers of VSe$_2$}
Moreover, we examine the energetics of the structural phases transition in the top layer of bi- and trilayers VSe$_2$ with different stacking models (Fig. \ref{fig8}). The notation of the types of Bernal stacking is similar to graphite. These results also can be applied for VSe$_2$ non-covalently attached to substrates.
\begin{figure}[ht]
    \includegraphics[width=1\columnwidth]{Fig8.png}
    \caption{Schematic representation of the unit cells used for simulating VSe$_2$ trilayers characterized by different stacking models.}
    \label{fig8}
\end{figure}

Results of the calculations (Fig. \ref{fig9}) evidence similarity the case of few-layer VSe$_2$ with monolayer. The configuration of the H type corresponds to the structural ground state for all types of the stacking in few-layer case. The energy required for the transition from T to H phase is about 0.60 eV for AA- and AB- stacking in bilayer. In trilayer the most energetically favorable stacking orders are AAA and ABC.
\begin{figure}[ht]
    \centering
    \includegraphics[width=1\columnwidth]{Fig9.png}
    \caption{Total energy (left panels) and magnetic moment (right panels) of two- (a,b) and three-layer (c,d) VSe$_2$ systems estimated for H, T and intermediate structures.}
    \label{fig9}
\end{figure}

 Thus, similarly to free standing VSe$_2$ monolayer, in top layer of VSe$_2$ there can be transition between two structural configurations at a moderate heating. The magnetic moments of vanadium atoms belonging to the upper layer of the few-layer structures change from 0.64$\mu_B$ to 0.82$\mu_B$. Such a change is fully connected with a redistribution of the electrons between different spin channels, main contributions are from $xy$ and $x^2-y^2$ orbitals of vanadium atoms similar to monolayer case. Therefore, the presence of the substrates does not influence significantly on sensitivity of bi- and trilayer VSe$_2$ systems to structural changes.

\subsection{Interlayer binding energy}
To understand the effect of interlayer interactions on structural properties we have
checked interlayer distances and binding energies. The binding energies E$_b$ for
different VSe$_2$ structures were calculated by using the following expression E$_b=($E$-n*$E$_{mono})/m$, where E is the total energy of considered system, E$_{mono}$ - total energy of monolayer, $n$ - number of layers in the considered system, $m$ - average number of interlayer interactions ($m$ = 2, 3/2 and 1 for bulk, 3 and 2-layers,
respectively). Results of these calculations are presented in Table
\ref{tab1}.

In the case when van der Waals interaction is neglected we obtain that the distance between V-V atoms belonging to the same layer is 3.33 \AA \ and the Se-Se interlayer distance equals to 3.12 \AA. When the van der Waals interaction is taken into account such distances equal to 3.31 \AA \ and 3.04 \AA, respectively. In 2- and 3-layer cases we considered the structures (Fig.\ref{fig8}) with the lowest total energies. Calculated values of the binding energies evidence that few-layer VSe$_2$ is pure van der Waals structure in contrast to bulk VSe$_2$ where London dispersion forces is a small addition to electrostatic interactions between V-cations and Se-anions from different
layers. The changes of interlayer distances are proportional to contribution of the dispersion forces to the binding energies (about 0.1 \AA \, in bulk and 0.4 - 0.6 \AA \, in few-layer systems). Therefore, the energy difference in migration barriers in bulk and few-layer VSe$_2$
can be explained by contribution from electrostatic repulsion of anions from the layer
above.
\begin{footnotesize}
\begin{table}[h!]
\begin{tabular}{c c c c  }
\hline
VSe$_2$  & E$_b$  with & E$_b$  without & Interlayer distance \\
structure & vdW, meV&  vdW, meV &with vdW \\
 & & &(without vdW), \AA\\
\hline
 T-bulk& 7.93&4.79&3.04(3.12)\\
 H-bulk&99.67&95.78&3.22(3.32)\\
 T-two&15.51&-9.36&3.11(3.51)\\
 H-two&24.68&-29.86&3.69(4.27)\\
 T-three&19.05&-90.86&3.08(3.57)\\
 H-three&100.82&-51.77&3.66(4.14)\\
\hline
\end{tabular}
\caption{Interlayer binding energies (meV/formula unit) and interlayer distances calculated for different VSe$_2$ structures with and without vdW interaction.}
\label{tab1}
\end{table}
\end{footnotesize}

\subsection{Phonon dispersion}
To complete the picture of physical properties of VSe$_2$ monolayer we have performed calculations of phonon dispersions by using VASP and Phonopy packages \cite{phonopy}. These
combination of the packages is widely used for studying of vibrational properties in
similar systems.\cite{Ersan} For such calculations we used $3\times3\times1$ supercell to obtain sets of forces and mesh grids: $10\times10\times1$ for monolayer and $6\times6\times6$ for bulk. Both H and T phases in nonmagnetic and ferromagnetic configurations were considered.

\begin{figure}[ht]
    \centering
    \includegraphics[width=1\columnwidth]{Fig10.png}
    \caption{Phonon dispersions calculated for nonmagnetic  (red dashed line) and ferromagnetic state (blue solid line) of monolayer and bulk VSe$_2$. Both T and H phase structures are presented. Red dots denote experimental frequencies taken from Ref.\cite{phonon1}.}
    \label{fig10}
\end{figure}
The calculated phonon spectra are presented in Fig. \ref{fig10}. For the T phase systems (bulk and monolayer) the resulting dispersions demonstrate a weak sensitivity to the magnetism. It is not the case for the H phase configurations. In the nonmagnetic state for H-monolayer and H-bulk we observe a soft phonon mode in the direction $\Gamma-$M$-$K for monolayer and in all symmetry directions for bulk. Existence of such a mode indicates structural instability. Importantly, in the ferromagnetic case the soft mode disappears, which means that the account of magnetism provides structural stability of H phase in both monolayer and bulk. The cause of this effect of magnetic configurations is the robustness of
magnetic interactions (see Fig.\ref{fig6} and discussion above) which is
the same order of magnitude as difference between structural
phase.  For H-bulk and T-monolayer ferromagnetic systems the calculations reveal the appearance of the indirect gap of ~0.57 THz. Comparison of the calculated dispersion curves with available experimental data from Ref.\cite{phonon1} obtained by a point-contact spectroscopy and Raman methods can be fulfilled only for the $\Gamma$ point for which experimental oscillation frequencies are ~6.04 (25 meV) and ~9.67 (40 meV) THz. Our theoretical values of 6.28 and 10.42 THz are in good agreement with experimental data.

\subsection{Structural phase transition by stretching}
The last step of our survey is the modeling of stretch which can appear in the monolayer due to substrate influence. To simulate  this effect, we increase $a$ and $b$ lattice vectors of our structure and then relax atomic positions to find a new ground state corresponding to new lattice parameters. Results of the calculation evidence that a stretching more than 3 percent leads to phase transition of the ground state configuration from H to T in monolayer and bilayer VSe$_2$ (Fig. \ref{fig11}a).
\begin{figure}[ht]
    \centering
    \includegraphics[width=1\columnwidth]{Fig11.png}
    \caption{Energy difference between H and T structures of VSe$_2$ (a) and energy barrier (b) as functions of stretching in $a$ and $b$ lattice directions. Lines of different colors correspond to systems with different numbers of layers.}
    \label{fig11}
\end{figure}

Therefore, the experimentally observed structure\cite{xmcd} of T type can result from a substrate-induced strain. Another effect of the stretching is a decreasing energy barrier for migration between different configurations (Fig. \ref{fig11}b). Here we define the energy barrier as energy difference between T structure and intermediate 30$^{\circ}$ structure

\section{CONCLUSIONS}
Results of first-principles calculations demonstrate that the energy barrier for the transition between two structural states of VSe$_2$ monolayer with a step-by-step rotation of the single Se atom is rather high. From the other hand the energy cost of the rotation of whole selenium layer is rather low (about 0.60 eV for monolayer and 0.80 eV for bulk). In the case of the monolayer it could be realized with a heating of the samples. The excitation energies of the rotation of the selenium layer up to 10$^\circ$ are very low, therefore, the realistic theoretical description of VSe$_2$ (from monolayer to bulk) should take into account these small deviations from ideal crystal structure.

Our calculations demonstrate that the transition from the experimentally observed T configuration to the H configuration is accompanied by a considerable change in electronic structure which is a redistribution of $3d$ electrons of vanadium between orbitals. Such transitions significantly influence on transport and thermal properties of VSe$_2$. From the other hand, the values of magnetic moments and total energies of ferro- and antiferromagnetic configurations change gradually between two structural phases.

In all the considered cases (bulk, few-layer and monolayer) system demonstrate strong favorability of ferromagnetic structure. The analysis of the calculated phonon dispersions has demonstrated a principal role of the ferromagnetism in stabilization of the atomic structure of the VSe$_2$ monolayer in H phase and similar systems. On the basis of the obtained results we can conclude that the experimentally observed paramagnetism in bulk VSe$_2$ and contradictory results of magnetic measurements for monolayers on different substrates are not caused by structural changes.

The calculations for bi- and trilayers demonstrate that the energy barrier of transition is similar to monolayer. The strain, possibly induced by the substrate, provides the change of the most energetically favorable structure from H to T. Therefore, the experimental observation of T configuration can result from a VSe$_2$ structure stretching by more than 3 percent on substrates. Another effect of the stretching is a decrease of the energy barrier of transition between structural phases. Thus both strain and deviation from ideal structure should be taken into account for realistic description of VSe$_2$ monolayer on substrates.

\section*{Conflicts of interest}
There are no conflicts to declare.

\section*{Acknowledgements}
This work was supported by the Russian Science Foundation, Grant No. 18-12-00185.

\balance


\section{Introduction and Notation}

Principal component analysis (PCA) is a standard tool for dimensionality reduction. Motivated by its extensions to functional PCA and kernel PCA, we are concerned with the statistical properties of PCA in infinite dimensions. In such scenarios, the eigenvalues of the covariance (resp. kernel) operator often decay at a polynomial or nearly exponential rate, which in turn implies that the minimal reconstruction error has a certain rate in the reconstruction dimension. In this paper, we investigate whether the latter rate is also achieved by the empirically chosen model.

We consider a random variable $X$ with values in a Hilbert space $\mathcal{H}$. In what follows, $\langle\cdot,\cdot\rangle$ denotes the inner product in $\mathcal{H}$, with $\norm{\cdot}$ being the corresponding norm. We suppose that $X$ is  centered and strongly square integrable, meaning that $\mathbb{E}X=0$ and $\mathbb{E}\| X\| ^2< \infty$. The goal is to reduce the dimensionality of $X$ by finding the minimizer $P_{\leq d}$ of the reconstruction error $R(P)=\mathbb{E}\| X-PX\| ^2$ over the class $\mathcal{P}_d$ of orthogonal projections of rank $d$. Yet, the distribution of $X$ is unknown and therefore the minimizer $P_{\leq d}$ cannot be computed. The idea is, for a given sequence $X_1, \dots, X_n$ of independent copies of $X$, to compute the minimizer $\hat{P}_{\leq d}\in \mathcal{P}_d$ of the empirical reconstruction error $R_n(P)= n^{-1}\sum_{i=1}^n \| X_i - PX_i\| ^2$. Given the interest in the performance of the empirically chosen model $\hat{P}_{\leq d}$ especially when observing new data, it is natural to analyze the reconstruction error of PCA, i.e.~the random variable $R(\hat{P}_{\leq d})$.

Bounds for the reconstruction error of PCA and the corresponding excess risk $R(\hat{P}_{\leq d})-R(P_{\leq d})$ can be derived using the theory of empirical risk minimization. This has
been pursued by Shawe-Taylor, Williams, Cristianini and Kandola \cite{SWCK05} and Blanchard, Bousquet and Zwald \cite{BBZ07}. In \cite{SWCK05}, a slow $n^{-1/2}$-rate is derived, while
in \cite{BBZ07}, it is shown that the convergence rate of the excess risk can be faster than $n^{-1/2}$ (but depending on a spectral gap condition).
In Rei\ss{} and Wahl \cite{ReissWahl}, it is shown how a fast $n^{-1}$-rate can be obtained by a version of the well-known Davis--Kahan $\sin\Theta$ theorem.
Moreover, since the resulting bound breaks down quickly in the infinite-dimensional case, \cite{ReissWahl} and subsequently Jirak and Wahl \cite{MR4052188} and Wahl \cite{Wahl} developed new perturbation techniques, leading (among others) to sharp excess risk bounds in the case that the
eigenvalues of the covariance operator of $X$ have exponential or polynomial decay.

Many algorithms in statistics and machine learning employ PCA as a first step to reduce the high dimensionality of the data. Analyzing such procedures, it often
turns out that the main interest is not in optimal rates for the excess risk but rather in bounds implying that empirical PCA produces a finite-dimensional model
having (up to a constant) the same reconstruction property as the optimal one (see, for example, Nouy \cite{Nouy}). Furthermore, the dependence on PCA is often
highly non-linear, hence expectation bounds, as derived in \cite{ReissWahl}, are not sufficient to analyze these procedures. Motivated by these two facts, our
main goal is to derive oracle inequalities in high-probability that can be expressed relative to the minimal reconstruction error. In particular, a major
achievement of our paper is to show that oracle inequalities of the type
\begin{equation}\label{weakOracleBound}
R(\hat{P}_{\leq d})\leq C \min_{P \in \mathcal{P}_d} R(P)
\end{equation}
can be established (with high probability) under one-sided eigenvalue conditions. For instance, we apply our bounds to eigenvalues satisfying polynomial or exponential upper bounds,
which are key features in reproducing kernel Hilbert spaces and functional data approaches in machine learning and statistics, see
e.g.~\cite{BachJordan,MR2335249,pmlr-v75-belkin18a}.

We finish the introduction by describing the link of PCA to the spectral decomposition of the (empirical) covariance operator. The covariance operator of $X$ is denoted by $\Sigma=\mathbb{E}X\otimes X$. By the spectral theorem, there exists a sequence $\lambda_1\geq \lambda_2\geq\dots>0$
of  positive eigenvalues (which is summable since $\mathbb{E}\| X\| ^2< \infty$) together with an orthonormal system of eigenvectors $u_1,u_2,\dots$ such that $\Sigma=\sum_{j\ge 1}\lambda_j P_j$, with rank-one projectors $P_j=u_j\otimes u_j$,
where $(u\otimes v)x=\langle v,x\rangle u$, for $u,v,x\in \mathcal{H}$. Without loss of generality we shall assume that the eigenvectors $u_1,u_2,\dots$ form an orthonormal basis of $\mathcal{H}$ such that $\sum_{j\ge 1}P_j=I$. In addition, we define the sample covariance of $X_1, \dots, X_n$ as $\hat{\Sigma} = n^{-1} \sum_{i = 1}^n X_i \otimes X_i$. Again, there exists a sequence $\hat{\lambda}_1 \geq \hat{\lambda}_2 \geq \cdots \geq 0$ of eigenvalues of $\hat{\Sigma}$ and an orthogonal basis of eigenvectors $\hat{u}_1, \hat{u}_2, \dots$ of $\mathcal{H}$ such that $\hat{\Sigma} = \sum_{j \geq 1} \hat{\lambda}_j \hat{P}_j$ with $\hat{P}_j = \hat{u}_j \otimes \hat{u}_j$. Now, the reconstruction error and the empirical reconstruction error can be written as
\[R(P) = \textup{tr}(\Sigma( I-P))\quad\text{and}\quad R_n(P) =\textup{tr}(\hat{\Sigma} (I-P)),
\] where $\textup{tr}(\cdot)$ denotes the trace. Hence, a minimizer of $R(P)$ maximizes $\textup{tr} (\Sigma P)$, from which one can deduce that a minimizer of the reconstruction error is given by the orthogonal projection onto the linear subspace spanned by the first $d$ eigenvectors of $\Sigma$, i.e.~\[
P_{\leq d} = \sum_{j \leq d} P_j\quad\text{satisfies}\quad P_{\leq d} \in \arg\min_{P \in \mathcal{P}_d}R(P)\quad\text{with}\quad R(P_{\leq d}) = \sum_{k>d} \lambda_k.
\]
Replacing the projections $P_j$ by $\hat{P}_j$, the minimizer $\hat{P}_{\leq d}$ of the empirical reconstruction error is given in the same way.

\section{Main result}

In this section we formulate our main error bound for the reconstruction error of PCA. It relies on a sub-Gaussian assumption on the Karhunen--Lo\`eve coefficients $\eta_j=\lambda_j^{-1/2}\langle X,u_j\rangle$, $j\geq 1$. Extensions under weaker moment assumptions are possible but beyond the scope of this paper.
\begin{assumption}\label{AssSG}
 Suppose that the $(\eta_j)_{j \geq 1}$ are independent. Moreover, for some constant $L > 0$ suppose that
\begin{equation}
 \sup_{j\geq 1}\sup_{q\geq 1} q^{-1/2} \left(\mathbb{E}|\eta_j| ^q\right)^{1/q} \leq L.
\end{equation}
\end{assumption}
Assumption \ref{AssSG} is satisfied with $L=1$ if $X$ is Gaussian. It gives a slightly stronger notion of a sub-Gaussian random variable than the
one in e.g. \cite[Definition 5.7]{Vershynin} or \cite[Assumption 2.1]{ReissWahl}.
\begin{theorem}\label{Theorem}
Under Assumption \ref{AssSG}, there are constants $c_1,c_2,C_1>0$ depending only on $L$ such that the following holds. Let $1\leq d'\leq d$ be natural numbers such that $\lambda_{d'}\geq 2\lambda_{d+1}$ and $\max(d',\lambda_{d'}^{-1}\sum_{k>d} \lambda_k) \leq c_1 n$. Then, for all $1\leq t \leq c_2n$, with probability at least $1-\exp(-t)$,
\begin{equation}\label{EqTheorem1}
 R(\hat{P}_{\leq d}) \leq \left( 1 + C_1 \left(\frac{d'}{n}+\frac{t}{n}\right)\right) \min_{P\in \mathcal{P}_{d'}}R(P).
\end{equation}
In particular, with probability at least $1-\exp(-c_2n)$,
\begin{equation}\label{EqTheorem2}
  R(\hat{P}_{\leq d}) \leq (1+C_1(c_1+c_2))\min_{P \in \mathcal{P}_{d'}} R(P).
\end{equation}
\end{theorem}
The second bound \eqref{EqTheorem2} shows that the reconstruction error of PCA can be bounded by a constant (typically close to one) times the minimal reconstruction error over $\mathcal{P}_{d'}$.

\section{Examples and discussion}
\label{SecEx}
Let us illustrate our upper bound for eigenvalues satisfying a polynomial or nearly exponential upper bound. Such eigenvalue structures are typically considered
in the context of functional data or statistical machine learning, see e.g. \cite{BachJordan,MR15, pmlr-v75-belkin18a}. In these cases, $\min_{P \in \mathcal{P}_d} R(P)$ (equals the
remainder trace) has a certain rate in $d$ and our bounds imply that the same rate holds for $R(\hat{P}_{\leq d})$.
\begin{corollary}
\label{CorPD}
Suppose that Assumption \ref{AssSG} holds. Moreover, suppose that for some $\alpha>1$ and $K>0$ we have $\lambda_j \leq K j^{-\alpha}$ for all $j\geq 1$. Then there are constants $c_1,c_2,C_1>0$ depending only on $\alpha,K$ and $L$ such that, for all $d\leq c_1 n$, with probability at least $1-\exp(-c_2n)$,
\begin{equation}\label{EqCorPD}
  R(\hat{P}_{\leq d}) \leq  C_1 d^{1-\alpha}.
\end{equation}
\end{corollary}
If additionally $\lambda_j \geq K^{-1} j^{-\alpha}$ for all $j\geq 1$, then we have $\min_{P\in \mathcal{P}_{d}}R(P)=\sum_{k>d}\lambda_k\geq c_3d^{1-\alpha}$ for a constant $c_3>0$, and \eqref{EqCorPD} can be reformulated as $R(\hat{P}_{\leq d}) \leq C_1\min_{P\in \mathcal{P}_{d}}R(P)$.
\begin{corollary}
\label{CorNED}
Suppose that Assumption \ref{AssSG} holds. Moreover, suppose that for $\alpha,K>0$ and $\beta\in (0,1]$ we have $K^{-1}\exp(-\alpha j^\beta)\leq \lambda_j \leq K \exp(-\alpha j^\beta)$ for all $j\geq 1$. Then there are constants $c_1,c_2,C_1,C_2>0$ depending only on $\alpha,\beta,K$ and $L$ such that, for all $d\leq c_1 n$, with probability at least $1-\exp(-c_2n)$,
\[
R(\hat{P}_{\leq d}) \leq  C_1 d^{1-\beta}\exp(-\alpha (d+1)^\beta)\leq C_2\min_{P \in \mathcal{P}_{d}} R(P).
\]
\end{corollary}
Let us compare these bounds to the Davis--Kahan $\sin \Theta$ theorem, one of the most well-known perturbation bounds for eigenspaces; see e.g. \cite{MR3371006} for a recent
statistical account. Applied to the excess risk (cf. Proposition 2.2 in \cite{ReissWahl}) it gives for every $1\leq t\leq n$, with probability at least
$1-\exp(-t)$,
\begin{equation}
\label{EqDK}
        R(\hat P_{\le d})- \min_{P\in \mathcal{P}_d} R(P) \leq  C\frac{t}{(\lambda_d-\lambda_{d+1})n},
\end{equation}
where $C >0$ is a constant depending only on $L$ and the trace of $\Sigma$. We observe that the upper bound typically increases for fixed sample size $n$ and increasing $d$, whereas the minimal reconstruction error decreases in $d$. For instance,
if $\lambda_j = K \exp(-\alpha j^\beta)$, $j\geq 1$, the bound \eqref{EqDK} breaks down for $d$ of size $(\log n)^{1/\beta}$. In contrast,  Corollary \ref{CorNED} does not depend on any gap condition and gives sharp bounds for $d\leq c_1n$.

\section{Proof of  Theorem \ref{Theorem}}
Throughout the proof, we use the letters $c,C>0$ for constants depending only on $L$ that may change from line to line. We start with formulating a more technical version of our main result in terms of the weighted covariance operator
\[
\Sigma_{d'}:= S_{\leq d'}\Sigma S_{\leq d'}\quad \text{where} \quad S_{\leq d'} = \sum_{j\leq d'} (\lambda_j - \lambda_{d+1})^{-1/2}P_j  \quad\text{and}\quad d'\leq d.
\]
Note that $S_{\leq d'}$ can be interpreted as the square-root of a (partial) reduced resolvent.
For this operator it holds that
\begin{align}\begin{aligned}
\| \Sigma_{d'}\| _{\infty} =& \frac{\lambda_{d'}}{\lambda_{d'} - \lambda_{d+1}},\quad\textup{tr}(\Sigma_{d'}) = \sum_{j \leq d'} \frac{\lambda_j}{\lambda_j - \lambda_{d+1}},\\
&\| \Sigma_{d'}\| ^2_{2} = \sum_{j \leq d'} \frac{\lambda_j^2}{(\lambda_j-\lambda_{d+1})^2},
\end{aligned}\end{align}
where $\| \cdot\| _{\infty}$ and $\| \cdot\| _{2}$ denote the operator norm and the Hilbert--Schmidt norm, respectively. Using these quantities, we can express our main result as follows.
\begin{theorem}
\label{techThm}
Under Assumption \ref{AssSG}, there are constants $c_1',c_2',C_1'>0$ depending only on $L$ such that the following holds. For all natural numbers $d' \leq d$ satisfying
\begin{equation}
 \label{AssEV}
\| \Sigma_{d'}\| _{\infty} \textup{tr}(\Sigma_{d'})+\| \Sigma_{d'}\| _{\infty}\sum_{k>d} \frac{\lambda_k}{\lambda_{d'}-\lambda_k} \leq c_1'n,
\end{equation}
we have, for all $1\leq t\leq c_2'n/\| \Sigma_{d'}\| _{\infty}^2$, with probability at least $1-\exp(-t)$,
\begin{equation}
\label{EqtechThm}
R(\hat{P}_{\leq d}) \leq \left( 1 + C_1'\left(\| \Sigma_{d'}\| _{\infty}\textup{tr}(\Sigma_{d'})\frac{1}{n}+\| \Sigma_{d'}\| ^2_{\infty} \frac{t}{n}\right)\right) \sum_{k > d'} \lambda_k.
\end{equation}
\end{theorem}
Let us first show how  Theorem \ref{techThm} implies  Theorem \ref{Theorem}. First, Condition \eqref{AssEV} is implied by the assumptions of Theorem \ref{Theorem}, provided that $c_1=c_1'/8$, as can be seen by using that $\| \Sigma_{d'}\| _{\infty}\leq 2$ and $\textup{tr}(\Sigma_{d'}) \leq 2d'$.
Since additionally $\sum_{k > d'}\lambda_k = \min_{P\in \mathcal{P}_{d'}}R(P)$, \eqref{EqTheorem1} follows from \eqref{EqtechThm}.
To prove Theorem \ref{techThm}, we will need three technical statements.
The following perturbation bound follows from the proofs of Proposition 3.5 and Theorem 2.12 in Rei\ss{} and Wahl \cite{ReissWahl} (see equations (3.9) and
(3.14) applied with $\mu = \lambda_{d+1}$ and $r = s = d'$).
\begin{lemma}
\label{lem:boundRW}
Suppose that the assumptions of Theorem \ref{techThm} hold. Then, on the joint event
\begin{equation}
\label{EventRW}
 \{\| S_{\leq d'}(\Sigma - \hat{\Sigma})S_{\leq d'}\| _{\infty} \leq 1/4\} \cap \{\hat{\lambda}_{d+1} - \lambda_{d+1} \leq (\lambda_{d'}- \lambda_{d+1})/2\},
\end{equation}
we have
\[
 R(\hat{P}_{\leq d}) \leq 16 \| \Sigma_{d'}\| _{\infty}\| S_{\leq d'}(\Sigma - \hat{\Sigma})P_{>d'}\| ^2_2 + \sum_{k>d'}\lambda_k.
\]
\end{lemma}
The event in \eqref{EventRW} has been analyzed in Lemma 3.9 and Corollary 3.12 in Rei\ss{} and Wahl \cite{ReissWahl}, yielding the following lemma.
\begin{lemma}
\label{lem:pboundW}
Under the assumptions of Theorem \ref{techThm}, we have
{\small \begin{align*}
  \ensuremath\mathbb{P}&\left(\{\| S_{\leq d'}(\Sigma - \hat{\Sigma})S_{\leq d'}\| _{\infty} \leq 1/4\} \cap \{\hat{\lambda}_{d+1} - \lambda_{d+1} \leq (\lambda_{d'}- \lambda_{d+1})/2\}\right)\\
   &\hspace{7cm} \geq 1 - 2 \exp\left(-c_2'n/\| \Sigma_{d'}\| ^2_{\infty}\right).
\end{align*}}
\end{lemma}
The next step is to prove a concentration inequality for $\| S_{\leq d'} (\hat\Sigma - \Sigma)P_{>d'}\| _2$. It can be derived from a Hilbert space version of the Hanson--Wright inequality derived in Adamczak, Lata\l a and Meller
\cite{ALM}. Alternatively, one can also multiply out the Hilbert--Schmidt norm, resulting in a polynomial chaos of degree
$4$ in sub-Gaussian random variables, and then apply Theorem 1.4 in \cite{AdamczakWolff}.
\begin{lemma}
\label{LemCI}
Under the assumptions of Theorem \ref{techThm}, we have, for all $t> 0$, with probability at least $1-2\exp(-t)$,
\begin{equation}
\label{EqLemCI}
 \| S_{\leq d'} (\hat\Sigma - \Sigma)P_{>d'}\| _2^2\leq  C \left( \textup{tr}(\Sigma_{d'})\frac{1}{n} + \| \Sigma_{d'}\| _{\infty} \frac{t}{n} + \| \Sigma_{d'}\| _{\infty} \frac{t^2}{n^2}\right)\sum_{k>d'}\lambda_k.
\end{equation}
\end{lemma}

\begin{proof}
The main idea of the proof is that we can rewrite the left-hand side in \eqref{EqLemCI} in terms of the squared Hilbert--Schmidt norm of a quadratic form in the Karhunen--Lo{\`e}ve coefficients $\eta_{ij}:= \lambda_j^{-1/2} \langle X_i, u_j \rangle$ with values being Hilbert-Schmidt operators. Indeed, we have
\[
S_{\leq d'} (\Sigma - \hat{\Sigma})P_{>d'} = \sum_{j\leq d'} \sum_{k>d'} \sum_{i=1}^n \biggl\{ \frac{1}{n}\sqrt{\frac{\lambda_j\lambda_k}{\lambda_j -\lambda_{d+1}}}u_j \otimes u_k \biggr\} \eta_{ij} \eta_{ik}.
\]
By Assumption \ref{AssSG}, the $\eta_{ij}$, $(i,j)\in \{1,\dots,n\}\times \mathbb{N}:= \tilde{\mathbb{N}}$, are independent, centered and sub-Gaussian random variables, meaning that we can apply Corollary 16
in \cite{ALM} (with the norm being the Hilbert--Schmidt norm) to the quadratic form
$\sum_{(i_1,j_1), (i_2, j_2) \in \tilde{\mathbb{N}}}a_{(i_1,j_1), (i_2,j_2)}\eta_{i_1j_1}\eta_{i_2j_2}$ with Hilbert--Schmidt operators
\begin{align*}
a_{(i_1,j_1), (i_2,j_2)} &= \frac{1}{2n} \left\{ \sqrt{\frac{\lambda_{j_1}\lambda_{j_2}}{\lambda_{j_1} -\lambda_{d+1}}} (u_{j_1} \otimes u_{j_2}) \mathbf{1}_{\{j_1\leq d',j_2>d'\}}\right.\\
&\left.\hspace{1.1cm} +  \sqrt{\frac{\lambda_{j_2}\lambda_{j_1}}{\lambda_{j_2} -\lambda_{d+1}}}(u_{j_2}\otimes u_{j_1})\mathbf{1}_{\{j_1>d',j_2\leq d'\}}\right\}\mathbf{1}_{\{i_1=i_2\}}.
\end{align*}
Note that, while  Corollary 16 in \cite{ALM} is formulated for finite index sets, it extends to our setting using e.g.~that the $a_{(i_1,j_1),(i_2,j_2)}$ are summable.
Therefore, we consider $x = (x_{(i,j)})_{(i,j) \in  \tilde{\mathbb{N}}}$ such that each component $x_{(i,j)} \in \mathbb{R}$. Similarly, we define $(x_{(i_1, j_1), (i_2, j_2)})_{(i_1, j_1), (i_2, j_2) \in \tilde{\mathbb{N}}}$.\newpage Exploiting the particular diagonal structure of
$(a_{(i_1,j_1),(i_2,j_2)})$, we compute, using the Cauchy--Schwarz inequality,
\begin{align*}
U_1 :=& \sup_{\| x\| _2\leq 1}\sqrt{\sum_{(i_1,j_1)}\Big\|\sum_{(i_2,j_2)}a_{(i_1,j_1),(i_2, j_2)} x_{(i_2,j_2)}\Big\|_2^2}\\
=& \, \frac{1}{2n}\sqrt{\max\Big\{\textup{tr}(\Sigma_{d'})\lambda_{d'+1}\, , \, \| \Sigma_{d'}\| _{\infty} \sum_{k>d'}\lambda_k\Big\}},\\
U_2 :=& \sup_{\big\| (x_{(i_1, j_1), (i_2, j_2)})\big\|_2\leq 1} \Big\| \sum_{(i_1,j_1), (i_2,j_2)} a_{(i_1,j_1), (i_2,j_2)} x_{(i_1,j_1), (i_2,j_2)}\Big\|_2\\
=&\, \frac{1}{\sqrt{2n}}\sqrt{\| \Sigma_{d'}\| _{\infty}\lambda_{d'+1}},\\
V :=& \sup_{\| x\| _2\leq 1, \| y\| _2\leq 1} {\Big\| \sum_{(i_1,j_1), (i_2,j_2)} a_{(i_1,j_1), (i_2,j_2)} x_{(i_1,j_1)}y_{(i_2,j_2)}\Big\|}_2\\
 =&\, \frac{1}{\sqrt{2}}\frac{1}{n}\sqrt{\| \Sigma_{d'}\| _{\infty} \lambda_{d'+1}},
\end{align*}
where $\| x\|_2^2 = \sum_{(i,j)\in \tilde{\mathbb{N}}} x_{(i,j)}^2$ and $\| (x_{(i_1, j_1), (i_2, j_2)})\| _2^2 = \sum_{(i_1, j_1), (i_2, j_2) \in \tilde{\mathbb{N}}} x_{(i_1, j_1), (i_2, j_2)}^2$. Applying Corollary 16 from \cite{ALM} with $U:= U_1+U_2$ and $t > 0$, we get
\begin{align}\label{HansonWright}
\begin{aligned}
\ensuremath\mathbb{P} &\biggl( \| S_{\leq d'}(\Sigma- \hat{\Sigma})P_{>d'}\| _2 \geq t+C\sqrt{n^{-1}\textup{tr}(\Sigma_{d'})\sum_{k>d'}\lambda_k}\biggr)\\
&\hspace{5cm} \leq 2 \exp\left(-\frac{1}{C}\min\left\{\frac{t^2}{U^2}, \frac{t}{V}\right\}\right)
\end{aligned}
\end{align}
and thus, with probability at least $1-2\exp(-t)$,
\begin{align*}
&\| S_{\leq d'} (\Sigma- \hat{\Sigma})P_{>d'}\| _2^2\\
&\qquad \leq C \biggl(\sqrt{n^{-1}\textup{tr}(\Sigma_{d'})\sum_{k>d'}\lambda_k}+ \max\left\{U\sqrt{t}\, ,\, Vt\right\}\biggr)^2\\
&\qquad \leq C\left(\textup{tr}(\Sigma_{d'})\frac{1}{n}+\textup{tr}(\Sigma_{d'}) \frac{t}{n^2} + \| \Sigma_{d'}\| _{\infty} \frac{t}{n}+ \| \Sigma_{d'}\| _{\infty}\frac{t^2}{n^2}\right) \sum_{k>d'}\lambda_k.
\end{align*}
The second term can be dropped using \eqref{AssEV} yielding the claim.
\end{proof}

\begin{proof}[End of proof of Theorem \ref{techThm}]
We choose $0<t \leq c_2'\| \Sigma_{d'}\| _{\infty}^{-2} n$. Then $\| \Sigma_{d'}\| _{\infty} (t/n)^2$ in \eqref{EqLemCI} is dominated by $\| \Sigma_{d'}\| _{\infty} t/n$. Combining the results from Lemma \ref{lem:boundRW}--\ref{LemCI}, we can finally argue that \eqref{EqtechThm}
holds with probability of $1-4\exp(-t)$. Restricting  ourselves to $t\geq 1$, this inequality also holds with probability at least $1-\exp(-t)$ for all $1\leq t \leq c_2'\| \Sigma_{d'}\| _{\infty}^{-2}n$, provided that we adjust the constants $C_1'$ and $c_2'$ appropriately.
\end{proof}

\section{Proofs of the corollaries}
\label{AppProof}

\begin{proof}[Proof of Corollary \ref{CorPD}]
First, there is a constant $C>0$ depending only on $\alpha$ and $K$ such that $\sum_{k>d'}\lambda_k\leq C(d')^{1-\alpha}$ for every $d'\geq 1$. In what follows we let $c_1$ be the constant from Theorem \ref{Theorem} such that $d\leq c_1n$ and we assume that $C\geq 2K$. We now first consider
the case that there is a $d'\leq d$ such that $\lambda_{d'}\geq Cd^{-\alpha}$. Then we have $\lambda_{d'}\geq 2\lambda_{d+1}$, as well as $\max(d',\lambda_{d'}^{-1}\sum_{k>d} \lambda_k)\leq d\leq c_1n$. Hence, we can apply Theorem \ref{Theorem} to the maximal $d'\leq d$ with this property,
leading to $R(\hat{P}_{\leq d})\leq C_2\sum_{k>d'}\lambda_k< C_2Cd^{1-\alpha}+C_2\sum_{k>d}\lambda_k\leq 2C_2Cd^{1-\alpha}$ with probability at least $1-\exp(-c_2n)$ and some constant $C_2>0$ depending only on $L$. On the other hand, if $\lambda_{d'}<Cd^{-\alpha}$ for all $d'\leq d$
(implying that $\lambda_1<Cd^{-\alpha}$), then the trivial bound $R(\hat{P}_{\leq d})\leq \textup{tr}(\Sigma)$ leads to $R(\hat{P}_{\leq d})\leq Cd^{1-\alpha}+\sum_{k>d}\lambda_k\leq 2Cd^{1-\alpha}$. This completes the proof.
\end{proof}

\begin{proof}[Proof of Corollary \ref{CorNED}]
For all $d\geq 1$ we have
\begin{equation}
\label{EqCompRT}
C^{-1}(d+1)^{1-\beta}\exp(-\alpha (d+1)^\beta)\leq\sum_{k>d}\lambda_k\leq C d^{1-\beta}\exp(-\alpha d^\beta)
\end{equation}
with a constant $C>1$ depending only on $\alpha,\beta$ and $K$. In fact, this can be seen by a comparison of the sum with an integral combined with standard estimates for the incomplete Gamma function.  Moreover, by a concavity argument, we have
$\beta x\leq 1-(1-x)^\beta\leq 2^{1-\beta}\beta x$, $x\in[0,1/2]$ and $\beta\in(0,1]$, from which we deduce that for $d'=d+1-k$ and $k\leq (d+1)/2$, the inequality
\begin{equation}
\label{EqCompEV}
K^{-2}\exp(\alpha\beta (d+1)^{\beta-1}k)\leq \frac{\lambda_{d'}}{\lambda_{d+1}}\leqK^2\exp(2^{1-\beta}\alpha\beta (d+1)^{\beta-1}k)
\end{equation}
holds. We now verify the conditions of Theorem \ref{Theorem} with $d'=d+1-k$, $k=\lceil (\alpha\beta)^{-1}(d+1)^{1-\beta}\log(2K^2)\rceil$. In the following, we will assume without loss of generality that $k\leq (d+1)/2$ (meaning that $d$ is larger than a constant depending on $\alpha,\beta$ and
$K$), because the bound \eqref{EqTheorem1} is obvious in the opposite case. First, using \eqref{EqCompEV}, $k$ is chosen such that $\lambda_{d'}\geq 2\lambda_{d+1}$ holds. In addition, using \eqref{EqCompRT}, we have $\lambda_{d}^{-1}\sum_{k>d} \lambda_k \leq Cd^{1-\beta}$. Hence, the assumptions of Theorem
\ref{Theorem} are satisfied for $d\leq (c_1/C)n$, and \eqref{EqTheorem2} yields $ R(\hat{P}_{\leq d}) \leq C_2 \sum_{k>d'}\lambda_k$ with probability at least $1-\exp(-c_2n)$. Thus the first claim follows from inserting  \eqref{EqCompRT} for $d=d'$ followed by an application of \eqref{EqCompEV}.
Moreover, the second claim follows from the first one, by applying \eqref{EqCompRT}.
\end{proof}

\subsection*{Acknowledgement}
The research of Martin Wahl has been partially funded by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1294 ``Data Assimilation'', Project (A4) ``Nonlinear statistical inverse problems with random observations''.

\bibliographystyle{plain}