\section{Introduction}
Data-to-text (D2T) generation is the task of transforming structured data into a natural language text which represents it \citep{reiter2000building,gatt2018survey}.
The output text can be generated in several steps following a pipeline, or in an end-to-end (E2E) fashion. Neural-based E2E architectures recently gained attention due to their potential to reduce the human input needed for building D2T systems. A disadvantage of E2E architectures is the lack of intermediate steps, which makes it hard to control the semantic fidelity of the output \citep{moryossef2019astep,ferreira2019neural}.

We focus on a D2T setup where the input data is a set of RDF triples in the form of \textit{(subject, predicate, object)} and the output text represents \textit{all} and \textit{only} facts in the data. This setup can be used by all D2T applications where the data describe relationships between entities  \cite[e.g.][]{gardent2017webnlg,budzianowski_multiwoz_2018}.\footnote{The setup can be preceded by the \textit{content selection} for selecting the relevant subset of data \cite[cf.][]{wiseman_challenges_2017}.} In order to combine the benefits of pipeline and E2E architectures, we propose to use the neural models with a limited scope. We take advantage of three facts: (1) each triple can be lexicalized using a trivial template, (2) stacking the lexicalizations one after another tends to produce an unnatural sounding but semantically accurate output, and (3) the neural model can be used for combining the lexicalizations to improve the output fluency.

In traditional pipeline-based NLG systems \citep{reiter2000building}, combining the lexicalizations is a non-trivial multi-stage process. Text structuring and sentence aggregation are first used to determine the order of facts and their assignment to sentences, followed by referring expression generation and linguistic realization. We argue that with a neural model, combining the lexicalizations can be simplified as several iterations of \textit{sentence fusion}---a task of combining sentences into a coherent text \citep{barzilay2005sentence}.

Our contributions are the following:
\begin{enumerate}[label={\arabic*)},nosep]
\item We show how to reframe D2T generation as iterative text editing, which makes it independent of dataset-specific input data format and allows to control the output over a series of intermediate steps.
\item We perform initial experiments using our approach on two major D2T datasets (WebNLG and Cleaned E2E) and include a quantitative and qualitative analysis of the results.
\item We perform zero-shot domain adaptation experiments and show that our approach exhibits a domain-independent behavior.

\end{enumerate}

\begin{figure*}[t]
    \centering
    \includegraphics[trim={0.1cm 0.3cm 0.2cm 0},clip,width=\textwidth]{d2t_text_editing}
    \caption{An example of a single iteration of our algorithm for D2T generation. In Step 1, the template for the triple is selected and filled. In Step 2, the sentence is fused with the template. In Step 3, the result for the next iteration is selected from the beam by filtering and language model scoring.}\label{fig:alg}
\end{figure*}

\section{Background}

Improving the accuracy of neural D2T approaches has attracted a lot of research interest lately.
Similarly to us, other systems use a generate-then-rerank approach \citep{dusek_sequence--sequence_2016,juraska_deep_2018} or a classifier to filter incorrect output \citep{harkous2020have}.
\citet{moryossef2019bimproving,moryossef2019astep}
split the D2T process into a symbolic text-planning stage and a neural generation stage. Other works improve the robustness of the neural model \citep{tian2019sticking,kedzie_good_2019} or employ a natural language understanding model \citep{nie_simple_2019} to improve the faithfulness of the output. Recently, \citet{chen-etal-2020-shot} finetuned GPT-2 \citep{radford2019language} for a few-shot domain adaptation.

Several models were recently applied to generic text editing tasks. \textsc{LaserTagger} \citep{malmi2019encode}, which we use in our approach, is a sequence tagging model based on the Transformer \citep{vaswani2017attention} architecture with the BERT \citep{devlin-etal-2019-bert} pre-trained language model as the encoder. Other recent text-editing models without a pre-trained backbone include EditNTS \citep{dong-etal-2019-editnts} and Levenshtein Transformer \citep{gu2019levenshtein}.

Concurrently with our work, \citet{kale2020few} explored using templates for dialogue response generation.
They use the sequence-to-sequence T5 model \citep{raffel2019exploring} to generate the output text from scratch instead of iteratively editing the intermediate outputs, which leaves less control over the model.

\section{Our Approach}
We start from single-triple templates and iteratively fuse them into the resulting text while filtering and reranking the results. We first detail the main components of our system and then give an overall description of the decoding algorithm.

\subsection{Template Extraction}
We collect a set of templates for each predicate. The templates can be either handcrafted, or automatically extracted from the lexicalizations of the single-triple examples in the training data. For unseen predicates, we add a single fallback template: \textit{The <predicate> of <subject> is <object>.}

\subsection{Sentence Fusion}
We train an in-domain sentence fusion model.
We select pairs $(X, X')$ of examples from the training data consisting of $(n, n+1)$ triples and having $n$ triples in common. This leaves us with an extra triple $t$ present only in $X'$. To construct the training data, we use the concatenated sequence $X~\mathrm{lex}(t)$ as a source and the sequence $X'$ as a target, where $\mathrm{lex}(t)$ denotes lexicalizing the triple $t$ using an appropriate template.
As a result, the model learns to integrate $X$ and $t$ into a single coherent expression.

We base our sentence fusion model on \textsc{LaserTagger} \citep{malmi2019encode}. \textsc{LaserTagger}
is a sequence generation model which generates outputs by tagging inputs with edit operations: \texttt{KEEP} a token, \texttt{DELETE} a token, and \texttt{ADD} a phrase before the token. In tasks where the output highly overlaps with the input, such as sentence fusion, \textsc{LaserTagger} is able to achieve performance comparable to state-of-the-art models with faster inference times and less training data.

An important feature of \textsc{LaserTagger} is the limited size of its vocabulary, which consists of $l$ most frequent (possibly multi-token) phrases used to transform inputs to outputs in the training data. After the vocabulary is precomputed, all infeasible examples in the training data are filtered out. At the cost of limiting the number of training examples, this filtering makes the training data cleaner by removing outliers. The limited vocabulary also makes the model less prone to common neural model errors such as hallucination, which allows us to control the semantic accuracy to a great extent using only simple heuristics and language model rescoring.

\subsection{LM Scoring}
We use an additional component for calculating an indirect measure of the text fluency. We refer to the component as the \textsc{LMScorer}. In our case,  \textsc{LMScorer} is a pre-trained GPT-2 language model \citep{radford2019language} from the Transformers repository\footnote{\url{https://github.com/huggingface/transformers}} \citep{Wolf2019HuggingFacesTS} wrapped in the \textit{lm-scorer}\footnote{\url{https://github.com/simonepri/lm-scorer}} package. We use \textsc{LMScorer} to compute the score of the input text $X$ composed of tokens $x_1\ldots x_n$ as a geometric mean of the token conditional probability:
\begin{align*}
    \operatorname{score}(X) = \Bigg( \prod_{i=1}^{n}{P(x_i|x_1 \ldots x_{i-1})} \Bigg)^{\frac{1}{n}}.
\end{align*}

\begin{table}[t
    \centering\small
    \begin{tabular}{ll}
        \multicolumn{2}{c}{\textbf{WebNLG}}  \\[0.1cm]
     \textit{foundedBy} &         \pbox{10cm}{<obj> was the founder of <subj>.\\<subj> was founded by <obj>.} \\[0.2cm] \midrule  \\[-0.3cm]

        \multicolumn{2}{c}{\textbf{E2E (extracted)}} \\[0.1cm]
        \textit{area+food} &   \pbox{5.5cm}{<subj> offers <obj2> cuisine in the <obj1>. \\ <subj> in <obj1> serves <obj2> food.} \\[0.2cm] \midrule \\[-0.3cm]
        \multicolumn{2}{c}{\textbf{E2E (custom)}} \\[0.1cm]
        \textit{near} &   \pbox{5cm}{<subj> is located near <obj>. \\ <obj> is close to <subj>.}
    \end{tabular}
    \caption{Examples of templates we used in our experiments. The templates for the single predicates in the WebNLG dataset and the pairs of predicates in the E2E dataset are extracted automatically from the training data; the templates for the single predicates in E2E are created manually.}
    \label{tab:templates_ex}

\end{table}

\subsection{Decoding Algorithm}
\label{sec:alg}
The input of the algorithm (Figure~\ref{fig:alg}) is a set of $n$ ordered triples. First, we lexicalize the triple $t_0$ to get the base text $X_0$. We choose the lexicalization for the triple as the filled template with the best score from \textsc{LMScorer}. This promotes templates which sound more natural for particular values.
In the following steps $i=1\ldots n$, we lexicalize the triple $t_i$ and append it after $X_{i-1}$. We feed the joined text into the sentence fusion model and produce a beam with fusion hypotheses. We use a simple heuristic (string matching) to filter out hypotheses in the beam missing any entity from the input data. Finally, we rescore the remaining hypotheses in the beam with \textsc{LMScorer} and let the hypothesis with the best score be the base text $X_{i}$. In case there are no sentences left in the beam after the filtering step, we let $X_{i}$ be the text in which the lexicalized $t_i$ is appended after $X_{i-1}$ without fusion (preferring accuracy to fluency). The output of the algorithm is the base text $X_n$ from the final step.

\begin{table*}[t
    \centering\small
    \begin{tabular}{lccccc<{\hspace{2mm}}c>{\hspace{2mm}}ccccc} \toprule
     & \multicolumn{5}{c}{\bf WebNLG} & & \multicolumn{5}{c}{\bf Cleaned E2E} \\
    \cmidrule{2-6} \cmidrule{8-12}
     & {\it BLEU} & {\it NIST} & \hspace{-1mm}{\it METEOR}\hspace{-1mm} & \hspace{-1mm}{\it ROUGE$_L$}\hspace{-1mm} & \hspace{-1mm}{\it CIDEr}\hspace{-1mm} & & {\it BLEU} & {\it NIST} & \hspace{-1mm}{\it METEOR}\hspace{-1mm} & \hspace{-1mm}{\it ROUGE$_L$}\hspace{-1mm} & \hspace{-1mm}{\it CIDEr}\hspace{-1mm} \\
    {\bf baseline} & 0.277 & 6.328 & 0.379 & 0.524 & 1.614 & & 0.207 & 3.679 & 0.334 & 0.401 & 0.365 \\
    {\bf zero-shot } & 0.288 & 6.677 & 0.385 & 0.530 & 1.751 & & 0.220 & 3.941 & 0.340 & 0.408 & 0.473 \\
    {\bf w/fusion } & 0.353 & 7.923 & 0.386 & 0.555 & 2.515 & & 0.252 & 4.460 & 0.338 & 0.436 & 0.944 \\
    {\bf SFC } & 0.524 & - & 0.424 & 0.660 & 3.700 & & 0.436 & - & 0.390 & 0.575 & 2.000 \\
    {\bf T5 } & 0.571 & - & 0.440 & - & - & & -  & - & - & - & - \\ \bottomrule
    \end{tabular}
    \caption{Results of automatic metrics on the WebNLG and Cleaned E2E test sets. The comparison is made with the results from the papers on the Semantic Fidelity Classifier (SFC; \citealp{harkous2020have}) and the finetuned T5 model (T5; \citealp{kale2020text}).}
    \label{tab:results}

    \end{table*}

\section{Experiments}

\subsection{Datasets} The WebNLG dataset \citep{gardent2017webnlg} consists of sets of DBPedia RDF triples and their lexicalizations. Following previous work, we use version 1.4 from \citet{ferreira2018enriching}. The E2E dataset \citep{novikova-etal-2017-e2e} contains restaurant descriptions based on sets of attributes (slots). In this work, we refer to the cleaned version of the E2E dataset \citep{dusek_semantic_2019}. For the domain adaptation experiments, we use \textsc{DiscoFuse} \citep{geva-etal-2019-discofuse}, which is a large-scale dataset for sentence fusion.

\subsection{Data Preprocessing} For WebNLG, we extract the initial templates from the training data from examples containing only a \textit{single} triple. In the E2E dataset, there are no such examples; therefore our solution is twofold: we extract the templates for \textit{pairs} of predicates, using them as a starting point for the algorithm in order to leverage the lexical variability in the data (manually filtering out the templates with semantic noise), and we also create a small set of templates for each \textit{single} predicate manually, using them in the subsequent steps of the algorithm (this is possible due to the low variability of the predicates in the dataset).\footnote{In the E2E dataset, the data is in the form of key-value slots. We transform the data to RDF triples by using the name of the restaurant as a \textit{subject} and the rest of the slots as \textit{predicate} and \textit{object}. This creates \textit{n-1} triples for \textit{n} slots.} See Table~\ref{tab:templates_ex} for examples of templates we used in our experiments.

\subsection{Setup}
As a \emph{baseline}, we generate the best templates according to \textsc{LMScorer} without applying the sentence fusion (i.e.\ always using the fallback).

For the \emph{sentence fusion} experiments, we use \textsc{LaserTagger} with the autoregressive decoder with a beam of size 10. We use all reference lexicalizations and the vocabulary size $V=100$, following our preliminary experiments, which showed that filtering the references only by limiting the vocabulary size brings the best results (see Supplementary for details). We finetune the model for 10,000 updates with batch size 32 and learning rate $2 \times 10^{-5}$.
For the beam filtering heuristic, we check for the presence of entities by simple string matching in WebNLG; for the E2E dataset, we use a set of regular expressions from TGen\footnote{\url{https://github.com/UFAL-DSG/tgen}} \citep{dusek_semantic_2019}. We do not use any pre-ordering steps for the triples and process them in the default order.

Additionally, we conduct a \textit{zero-shot domain adaptation} experiment. We train the sentence fusion model with the same setup, but instead of the in-domain datasets, we use a subset of the balanced-Wikipedia portion of the \textsc{DiscoFuse} dataset. In particular, we use the discourse types which frequently occur in our datasets, filtering the discourse types which are not relevant for our use-case. See Supplementary for the full listing of the selected types.

\begin{table*}[ht] \small
    \begin{tabular}{l p{13.5cm}}
        \textbf{Triples}   & \textit{(Albert Jennings Fountain, deathPlace, New Mexico Territory); (Albert Jennings Fountain, birthPlace, New York City); (Albert Jennings Fountain, birthPlace, Staten Island)} \\ \midrule
       \textbf{Step \#0} & Albert Jennings Fountain died in New Mexico Territory. \\
       \textbf{Step \#1} & Albert Jennings Fountain, who died in New Mexico Territory, was born in \greenund{New York City}. \\
       \textbf{Step \#2} & Albert Jennings Fountain, who died in New Mexico Territory, was born in New York City, \greenund{Staten Island}. \\ \midrule
       \textbf{Reference} & Albert Jennings Fountain was born in Staten Island, New York City and died in the New Mexico Territory.
       \end{tabular}
    \caption{An example of correct behavior of the algorithm on the WebNLG dataset. Newly added entities are underlined, the output from Step \#2 is the output text.}\label{tab:ex0}
    \end{table*}

\section{Analysis of Results}

We compute the metrics used in the evaluation of the E2E Challenge \citep{duvsek2020evaluating}: BLEU \citep{papineni-etal-2002-bleu},  NIST \citep{doddington2002automatic}, METEOR \citep{banerjee-lavie-2005-meteor}, ROUGE$_L$ \citep{lin-2004-rouge} and CIDEr \citep{vedantam2015cider}. The results are shown in Table~\ref{tab:results}.
The scores from the automatic metrics lag behind the state-of-the-art, although both the fusion and the zero-shot approaches show improvements over the baseline. We examine the details in the following paragraphs, discussing the behavior of our approach, and we outline plans for improving the results in Section~\ref{sec:future}.

\paragraph{Accuracy vs.\ Variability} Our approach ensures zero entity errors, since the entities are filled verbatim into the templates and in case an entity is missing in the whole beam, a fallback is used instead. Semantic inconsistencies still occur, e.g.\ if a verb or function words are missing.

The fused sentences in the E2E dataset, where all the objects are related to a single subject, often lean towards compact forms, e.g.: \textit{Aromi is a family friendly chinese coffee shop with a low customer rating in riverside}. On the contrary, the sentence structure in WebNLG mostly follows the structure from the templates and the model performs minimal changes to fuse the sentences together. See Table~\ref{tab:ex0} and Supplementary for examples of the system outputs. Out of all steps, 28\% are fallbacks (no fusion is performed) in WebNLG and 54\% in the E2E dataset.
The higher number of fallbacks in the E2E dataset can be explained by a higher lexical variability of the references, together with a higher number of data items per example in the E2E dataset, making it harder for the model to maintain the text coherency over multiple steps.

\paragraph{Templates} On average, there are 12.4 templates per predicate in WebNLG and 8.3 in the E2E dataset. In cases where the set of templates is more diverse, e.g.\ if the template for the predicate \textit{country} has to be selected from \{\textit{<subject> is situated within <object>, <subject> is a dish found in <object>}\}, \textsc{LMScorer} helps to select the semantically accurate template for the specific entities. The literal copying of entities can be too rigid in some cases, e.g.\ \textit{Atatürk Monument (İzmir) is made of ``Bronze''}, but these disfluencies can be improved in the fusion step.

\paragraph{Reordering} \textsc{LaserTagger} does not allow arbitrary reordering of words in the sentence, which can limit the expressiveness of the sentence fusion model. Consider the example in Figure~\ref{fig:alg}: in order to create a sentence \textit{English is spoken in Dublin, the capital of Ireland}, the model has to delete and re-insert at least one of the entities, e.g.\ \textit{English}, which has to be present in the vocabulary.

\paragraph{Domain Independence} The zero-shot model trained on \textsc{DiscoFuse} is able to correctly pronominalize or delete repeated entities and join the sentences with conjunctives, e.g.\ \textit{William Anders was born in British Hong Kong\underline{, and was} a member of the crew of Apollo 8}. While the model makes only a limited use of sentence fusion, it makes the output more fluent while keeping strong guarantees of the output accuracy.

\section{Future Work}
\label{sec:future}
Although the current version of our approach is not yet able to consistently produce sentences with a high degree of fluency, we believe that the approach provides a valuable starting point for controllable and domain-independent D2T generation. In this section, we outline possible directions for tackling the main drawbacks and improving the results of the model with further research.

Building a high-quality sentence fusion model, which lies at the core of our approach, remains a challenge \cite{lebanoff2020learning}. Our simple extractive approach relying on existing D2T datasets may not produce sufficient amount of clean data. On the other hand, the phenomena covered in the \textsc{DiscoFuse} dataset are too narrow for the fully general sentence fusion. We believe that training the sentence fusion model on a larger and more diverse sentence fusion dataset, built e.g. in an unsupervised fashion \cite{lebanoff-etal-2019-scoring}, is a way to improve the robustness of our approach.

Fluency of the output sentences may be also improved by allowing more flexibility for the order of entities, either by including an ordering step in the pipeline \citep{moryossef2019astep}, or by using a text-editing model that is capable of explicit reordering of words in the sentence \cite{mallinson2020felix}. Splitting the data into smaller batches (i.e. setting an upper bound for the number of sentences fused together) could also help to improve the consistency of results with a higher number of data items.

Our string matching heuristic is quite crude and may lead to a high number of fallbacks. Introducing a more precise heuristic, such as a semantic fidelity classifier \citep{harkous2020have}, or a model trained for natural language inference \citep{dusek2020nli} could help to promote lexical variability of the text.

Finally, we note that the text-editing paradigm allows to visualize the changes made by the model, introducing the option to accept or reject the changes at each step, and even build a set of custom rules on top of the individual edit operations based on the affected tokens. This flexibility could be useful for tweaking the model manually for a production system.

\section{Conclusions}
We proposed a simple and intuitive approach for D2T generation, splitting the process into two steps: lexicalization of data and improving the text fluency. A trivial lexicalization helps to promote fidelity and domain independence while delegating the subtle work with language to neural models allows to benefit from the power of general-domain pre-training. While a straightforward application of this approach on the WebNLG and E2E datasets does not produce state-of-the-art results in terms of automatic metrics, the results still show considerable improvements above the baseline. We provided insights into the behavior of the model, highlighted its potential benefits, and proposed the directions for further improvements.

\section*{Acknowledgements}
We would like to thank the anonymous reviewers for the relevant comments. The work was supported by the Charles University grant No.~140320, the SVV project No.~260575, and the Charles University project PRIMUS/19/SCI/10.

\bibliographystyle{acl_natbib}


\section{Introduction}

In musically creative tasks which rely on some prior analysis of a music signal, (e.g., musical audio beat tracking for remixing or beat-matching) the mean performance of an algorithm over large datasets is of less relevance than knowing its exact performance for a given piece. Indeed, within Creative MIR, high importance is placed on the extremely accurate analysis of music signals \cite{andersen2016ismir}. Errors in analysis can have an immediate negative impact on the musical result, as well as being a confounding factor when seeking to evaluate the core creative task. To mitigate this problem, we have two options: to exploit high-level user-input to guide how the analysis is undertaken \cite{sapinto2019cmmr} towards a very accurate analysis from the outset; and/or to manually correct the output of the algorithm using a tool such as Sonic Visualiser \cite{cannam2010acm} until it is satisfactory for the end-user. To reduce the workload and cognitive burden, we normally wish to minimize the number of user interactions required.

For musical audio beat tracking, the alteration of beat detections may be challenging due to the underlying difficulty of the musical material, but the correction process can be achieved using two simple editing operations: insertions and deletions -- combined with repeated listening to audible clicks mixed with the input. The number of insertions and deletions correspond to counts of \textit{false negatives} and \textit{false positives} respectively, and form part of the calculation of the well-known F-measure. While the F-measure is routinely used in beat tracking (and many other MIR tasks) to measure accuracy, we can also view it in terms of the
effort required to transform an initial set of beat detections to a final desired result (e.g., a ground truth annotation sequence). In this way, a high F-measure would imply low effort in manual correction and vice-versa.

In practice, correcting beat detections often relies on a third operation: the \textit{shifting} of poorly localised individual beats. This shifting operation is particularly relevant when correcting tapped beats, which can be subject to human motor noise as well as jitter and latency during acquisition. Under the logic of the F-measure calculation, shifting beat detections which fall outside tolerance windows are effectively counted twice: as a false positive \textit{and} a false negative. We argue that for beat tracking evaluation this creates a modest, but important, disconnect between common practice in annotation correction and a widely-used evaluation method. We propose that, wherever reasonable, the single operation of shifting should be prioritised over a deletion followed by an insertion.

\section{Measuring Annotation Efficiency}

We devise a straightforward calculation for the \textit{annotation efficiency} based on counting the number of shifts, insertions, and deletions. We provide an open-source Python implementation\footnote{\url{https://github.com/MR-T77/ShiftIfYouCan}} which graphically displays the minimum set and type of operations required to transform a sequence of initial beat detections in such a way as to maximize the F-measure when comparing the transformed detections against a ground truth annotation sequence. It is important to note that our goal is \textit{not} to transform the beat detections such that they are absolutely identical to the ground truth, (although such transformations are theoretically possible) but rather to perform as few operations as possible to ensure $F=1.00$, subject to a user-defined tolerance window.

\begin{figure*}[h!]
\centering
  \includegraphics[width=17cm]{Regular_vis.pdf}
  \includegraphics[width=17cm]{Double_vis.pdf}
  \caption{Visualisation of the operations required to transform beat detections to maximize the F-measure when compared to the ground truth annotations. (Top) Original beat detections vs. ground truth annotations. (Bottom) Double variation of beat detections vs. ground truth annotations. The inner tolerance window is overlaid on all annotations, where as the outer tolerance window is only shown for those detections to be shifted.}
  \label{fig:regular}
\end{figure*}

Around each ground truth annotation, we create an \textit{inner tolerance window} (set to $\pm70$\,ms) and count the number of \textit{true positives} (unique detections), $t^{+}$. We mark each matching detection and annotation pair as ``accounted for'' and remove them from further analysis. All remaining detections then become candidates for \textit{shifting} or \textit{deletion}. For each remaining annotation, we then look for the closest unaccounted for detection within an \textit{outer tolerance window} (set to $\pm1$\,s), which we use to reflect a localised working area for manual correction. If any detection exists, we mark it as a shift along with the required temporal correction offset. After the analysis of all unaccounted for annotations is complete, we count the number of shifts, $s$. Any remaining annotations correspond to false negatives, $f^{-}$, with leftover detections marked for deletion and counted as false positives, $f^{+}$. To give a measure of annotation efficiency we adapt the evaluation method in \cite{dixon2001jnmr} to include the shifts:
\begin{equation}
ae = t^{+} / (t^{+} + s + f^{+} + f^{-}).
\end{equation}
Reducing the inner tolerance window transforms true positives into shifts and thus sends $t$ and hence $ae$ to $0$. In this limit, the modified detections are then identical to the target sequence.

To allow for metrical ambiguity in beat tracking evaluation, it is common to create a set of variations of the ground truth by interpolation and sub-sampling operations. In our implementation, we flip this behaviour, and instead create variations of the detections. In this way, we can couple a global operation applied to all detections (e.g., interpolating all detections by two), with the subsequent set of local correction operations. Whichever variation has the highest annotation efficiency represents the shortest path to obtaining an output consistent with the annotations.

The fundamental difference of our approach compared to the standard F-measure is that we view the evaluation from a user workflow perspective, and essentially \textbf{we shift if we can}. By recording each individual operation, we can count them for evaluation purposes as well as visualising them, as shown in Figure.~\ref{fig:regular} which contrasts the use of the original beat detections compared to the double variation of the beats. The example shown is from the composition ``Evocaci\`{o}n'' by Jose Luis Merlin. It is a solo piece for classical guitar which features extensive \textit{rubato} and is among the more challenging pieces in the Hainsworth dataset \cite{Hainsworth2004}. By inspection we can see the original detections are much closer to the ground truth than the double variation. They require just two deletions (covering a pause in the performance), three shifts, and two insertions, and thus the annotation efficiency is much higher: $0.650$ compared to $0.361$.

The precise recording of the set of individual operations allows additional deeper evaluation which can indicate precisely which operations are most beneficial and in which order. For the F-measure, shifts are always more beneficial than the isolated insertions or deletions, but for other evaluation methods, i.e., those which measure continuity, the temporal location of the operation may be more critical. By viewing evaluation from a transformation perspective combined with an informative visualisation, we hope our implementation can contribute to a better qualitative understanding of beat tracking algorithms.

Looking beyond beat tracking, our approach can readily be applied to other temporal MIR tasks, e.g., onset detection and structural boundary detection -- subject to appropriate tolerance windows. In future work, we intend to expand the functionality to additionally correct labelled time instants, e.g., beat positions in each bar when estimating downbeats. We will also look to position our work in the context of the formal theory of edit distances. Finally, we plan to undertake a user-study to understand how well our measure of annotation efficiency correlates with user-reported effort in the manual correction of beat estimates.

\section{Acknowledgments}
This work is funded by national funds through the FCT - Foundation for Science and Technology, I.P., within the scope of the project CISUC - UID/CEC/00326/2020 and by European Social Fund, through the Regional Operational Program Centro 2020. It is also supported by the FCT-Foundation for Science and Technology, I.P., under the grant SFRH/BD/120383/2016 and the projects IF/01566/2015 and UIDP/00776/2020.


\section{Introduction} \label{Sec:Intro}
With the advance of technology, massive amounts of multivariate spatial data can be collected. As one example, researchers may use these datasets to investigate various issues in geographical, ecological (\cite{legendre2012numerical}) or atmospheric (\cite{von2001statistical}) sciences. Commonly, in these spatially correlated datasets there are dependencies both within and among the individual data processes, which makes modeling the multivariate spatial data a challenge. The difficulty of modeling is further intensified when the dimensionality $p$ is large. With a dataset of size $n$, it takes a total of $cp(p+1)/2$ parameters to describe the full covariance and cross-covariance structure of the model, where $c$ is the number of characteristic parameters per covariance and cross covariance. Further, it requires a computational cost of $\ensuremath{\mathcal{O}}(n^3p^3)$ for prediction using optimal linear predictors and for Gaussian likelihood evaluation
(\cite{cressie2015statistics,legendre2012numerical}).

One way to approach this problem is to use the spatial blind source separation (SBSS) framework see \cite{nordhausen2015blind} and \cite{bachoc2018spatial}. Blind source separation (BSS) is a well-studied multivariate procedure used to recover latent variables when only a linear mixture of them is observed; see for example \cite{comon2010handbook,NordhausenOja2018}.
 A common assumption for BSS is that the latent variables are second-order stationary and uncorrelated. That is, we assume $\bo x(\bo s)=\Om \z (\bo s)$, where $\bo x (\bo s) \in\R^p$ is the observed $p$-variate measurement at location $\bo s \in \mathbb{R}^d$, $\z (\bo s) \in\R^p$ is a latent second-order stationary $p$-variate source with uncorrelated components, and $\Om\in\R^{p\times p}$ is an unknown full-rank mixing matrix. To estimate the unmixing matrix $\G$, i.e. $\Om^{-1}$, \cite{nordhausen2015blind} have proposed an estimator based on the simultaneous diagonalization of two scatter matrices, and \cite{bachoc2018spatial} have extended this method to jointly diagonalize more than two scatter matrices for multivariate spatial data.

The SBSS model of \cite{nordhausen2015blind} gives no preference to any of the latent components, with all $p$ of them being basically of equal interest. However, in practical cases of BSS, it is often assumed that only a few components are of interest and to be regarded as the signal, while the remaining components are discarded as noise. That is, it is supposed that the latent components consist of two parts, $\z=(\z_{s}^T, \z_{w}^T)^T$, where $\z_{s}\in\R^q$ is the signal, and $\z_{w}\in\R^{p-q}$ is the noise. \cite{bootstrap_ts, virta2018determining} and \cite{NordhausenVirta2018} all consider components with serial dependence as signals in a time series context.

Following the work of \cite{virta2018determining}, in this study we consider SBSS in which signals are characterized  as components  having second-order spatial dependence.
We derive a test for the signal dimension $q$ based on the joint diagonalization of two or more scatter matrices that are specified by kernel functions, as done by \cite{bachoc2018spatial}. We then provide the asymptotic distribution of the test statistic. This asymptotic result does not rely on the strong assumption made by \cite{bachoc2018spatial}, that all the signal and noise components are asymptotically identifiable. To avoid this assumption, new proof techniques are needed, and we rely in particular on the extension of arguments made by \cite{virta2018determining} to a spatial setting. In addition, we demonstrate that using a different normalization factor for the scatter matrices than that used by \cite{bachoc2018spatial} enables the obtainment of a neater asymptotic distribution of the test statistic (see Remark \ref{remark:renormalization}). Based on the test, we then provide a consistent estimator of the unknown signal dimension.

We put forward several bootstrap versions of the test. For both the asymptotic and bootstrap tests, we demonstrate that computational gains are obtained when the observation locations are gridded. In an extensive simulation study, we then show that the various tests already have levels close to the nominal one, for small to moderate sample sizes. We also observe an accurate estimation of the signal dimension. We conclude that the asymptotic test is comparable to and often outperforms the bootstrap ones while being less computationally demanding. Employing an environmental application, we then show that our methods enable the reduction of the dimension of a multivariate spatial data set, retaining the most interpretable and informative estimated independent components and discarding the unusable ones as noise.

The remainder of the paper is organized as follows. In Section \ref{Sec:Setup}, we introduce the statistical setting of the problem and present our test statistic. The methods and main results are then described in Section \ref{Sec:Theory}, the simulation results are reported in Section \ref{Sec:Sim}, and the environmental application is outlined in Section \ref{Sec:Ex}. We finally discuss some concluding remarks in Section \ref{Sec:Conclusion}. The proofs of the theoretical results are provided in the appendix.

\section{Setup and Model}\label{Sec:Setup}

Suppose our data consists of a $p$-dimensional multivariate random field $\bo x(\bo s)=\{x_1(\bo s), \cdots, x_p(\bo s)\}^T$, $\bo s \in \mathcal{S}$, where $\mathcal{S}\subseteq\R^d$ is a region of interest. The covariance and cross-covariance functions of
$\bo x$, defining its second-order structure, are some of its central characteristics. We can refer, for instance, to \cite{deIaco2013using,genton2015cross} and \cite{gneiting2010matern} for an introduction and various approaches to modeling these.

Here, the second-order structure of $\bo x$ is assumed to obey an SBSS model:
\begin{equation}\label{eq:model}
\bo x(\bo s)=\Om \z(\bo s),
\end{equation}
where $\Om$ is a $p\times p$ unknown invertible matrix, and  $\z(\bo s)=\{z_1(\bo s),\cdots,z_p(\bo s)\}^T$ is the latent field having independent components with $\mathop{\mathrm{Cov}}\nolimits(\z(\bo s))=\I_p$ for all $\bo s\in\mathcal{S}$.

Let $I(\cdot)$ denote the indicator function throughout this paper and consider the kernel functions $f_0, f_1, \cdots, f_k$, with $f_\ell: \R^d\to\R$ for $\ell=0, \cdots, k$, and with $f_0(\bo s)=I(\bo s=\boldsymbol{0})$. For $f\in\{f_0, f_1, \cdots, f_k\}$, let
\[
F_{n,f}=\frac{1}{n}\sum_{i,j=1}^nf^2(\bo s_i-\bo s_j),
\]
where $\{ \bo s_1,\ldots,\bo s_n \} \subseteq \mathcal{S}$ is the set of two-by-two distinct observation points.
Notice that $F_{n,f_0}=1$. For $f\in\{f_1, \cdots, f_k\}$. The population local covariance (or scatter) matrices are then defined as,
\begin{equation}\label{eq:plc}
\bo M(f)=\frac{1}{n\sqrt{F_{n,f}}}\sum_{i=1}^n\sum_{j=1}^nf(\bo s_i-\bo s_j)\mathbb{E}\left(\bo x(\bo s_i)\bo x(\bo s_j)^T\right) \quad \text{and} \quad \bo M(f_0)=\frac{1}{n}\sum_{i=1}^n\mathbb{E}\left(\bo x(\bo s_i)\bo x(\bo s_i)^T\right),
\end{equation}
and the corresponding sample local covariance matrices are defined as
\begin{equation}\label{eq:slc}
\widehat{\bo M}(f)=\frac{1}{n\sqrt{F_{n,f}}}\sum_{i=1}^n\sum_{j=1}^nf(\bo s_i-\bo s_j)\bo x(\bo s_i)\bo x(\bo s_j)^T \quad \text{and} \quad \widehat{\bo M}(f_0)=\frac{1}{n}\sum_{i=1}^n\bo x(\bo s_i)\bo x(\bo s_i)^T.
\end{equation}
\begin{remark} \label{remark:renormalization}
The normalizing quantity $nF_{n,f}^{1/2}$ in \eqref{eq:plc} and \eqref{eq:slc} is slightly different from that in \cite{bachoc2018spatial}, where simply $n$ is used. Here, the introduction of $F_{n,f}^{1/2}$ enables us to obtain a simple and elegant asymptotic distribution of the test statistic for the number of noise components (see Proposition \ref{prop:asymt}).
\end{remark}

The $k+1$ sample local covariance matrices
$
\widehat{\bo M}(f_0), \widehat{\bo M}(f_1), \cdots, \widehat{\bo M}(f_k)
$
are used to estimate the unmixing matrix  $\G$ as
\begin{equation}\label{eq:hg}
\widehat{\bo \Gamma}\in\argmax_{\substack{\G: \G\widehat{\bo M}(f_0)\G^T=\I_p\\
\G\text{ has rows }\g_1^T,\cdots,\g_p^T\\
(\sum_{\ell=1}^k\{\g^T_j\widehat{\bo M}(f_{\ell})\g_j\}^2)_{j=1,\cdots,p}\text{ are in descending order}}} \sum_{\ell=1}^k\sum_{j=1}^p\{\g^T_j\widehat{\bo M}(f_{\ell})\g_j\}^2.
\end{equation}
The unmixing matrix should ``diagonalize" all $k$ local covariance matrices and we let for $\ell=1,\cdots,k$,
\[
\widehat{\bo D}_{\ell}=\widehat{\bo \Gamma}\widehat{\bo M}(f_{\ell})\widehat{\bo \Gamma}^T\]
where all $\widehat{\bo D}_{\ell}$ should be close to a diagonal matrix. Note that for finite data exact diagonalization is usually possible only for $k=1$. Further, by definition,
$
\sum_{\ell=1}^k\widehat{\bo D}^2_{\ell,1,1}\ge\cdots\ge\sum_{\ell=1}^k\widehat{\bo D}^2_{\ell,p,p}.$
We are now interested in the case in which there are $q$ ``real" continuous random fields in $\z$, while the remaining $p-q$ components are white noise.

For $q\in\{0,\cdots,p-1\}$, we are interested in testing the following hypothesis
\begin{center}
$H_{0q}:$ There are exactly $p-q$ white noise processes in $\z$.
\end{center}

This hypothesis is formalized in the following two conditions:
\begin{condition}\label{c:1}
For $a=1, \cdots, p-q$, the covariance function of $z_{q+a}$ is given by
$$\mathop{\mathrm{Cov}}\nolimits(z_{q+a}(\bo u), z_{q+a}(\bo v))=I(\bo u-\bo v=\boldsymbol{0}).$$
\end{condition}
\begin{condition}\label{c:2}
For $\ell=1,\cdots,k, f_{\ell}$ is symmetric and satisfies $f_{\ell}(\bm{0})=0$.
For $a=1, \cdots, q$, we have
\[
\liminf_{n\to\infty}\sum_{\ell=1}^k\bo M(f_{\ell})^2_{a,a}>0.\]
\end{condition}
Note that in Conditions \ref{c:1} and \ref{c:2}, we assume that the sources are ordered such that the $q$ signal components come first and are followed by the $p-q$ noise components. As the order of the sources is not identifiable, this assumption comes without loss of generality. The fulfillment of Condition \ref{c:2} means that the correlation in the signal fields $z_1,\cdots,z_q$ is sufficient for these signals to be asymptotically separated from the noise fields $z_{q+1}, \cdots, z_{p}$. It should also be noted that we do not need to consider the stronger assumption that the $q$ vectors of the $a$-th diagonal elements in $\bo M(f_1), \cdots, \bo M(f_k)$, for $a=1,\cdots,q$, for the signal random fields, are asymptotically distinct (see Assumption 9 in \cite{bachoc2018spatial}) and non-zero. Conditions \ref{c:1} and \ref{c:2} motivate the following block decompositions for $\ell=1,\cdots,k$:
\[
\widehat{\bo M}(f_{\ell})=\begin{pmatrix}
\widehat{\bo M}(f_{\ell})_{qq} & \widehat{\bo M}(f_{\ell})_{q0}\\
\widehat{\bo M}(f_{\ell})_{0q} & \widehat{\bo M}(f_{\ell})_{00}
\end{pmatrix} \quad \text{ and } \quad
\widehat{\bo D}_{\ell}=\begin{pmatrix}
\widehat{\bo D}_{\ell,qq} & \widehat{\bo D}_{\ell,q0}\\
\widehat{\bo D}_{\ell,0q} & \widehat{\bo D}_{\ell,00}
\end{pmatrix}
\]
where the blocks $\widehat{\bo M}(f_{\ell})_{qq}$ and $\widehat{\bo D}_{\ell,qq}$ have size $q\times q$ and the blocks $\widehat{\bo M}(f_{\ell})_{00}$ and $\widehat{\bo D}_{\ell,00}$ have dimension $(p-q)\times (p-q)$.\\

Then our test statistic is
\begin{equation}\label{eq:t}
t_q=\frac{n}{2}\sum_{\ell=1}^k||\widehat{\bo D}_{\ell,00}||^2,
\end{equation}
where $||\cdot||$ is the Frobenius norm. The test statistic is then expected to be bounded under the null hypothesis, and to diverge when one of $z_{q+1},\ldots,z_{p}$ is spatially correlated. The test will reject the null hypothesis $H_{0q}$ if $t_q$ is larger than a certain threshold, in which case the dataset provides indications that more than $q$ signal components are present.  For a nominal level $\alpha \in (0,1)$, the threshold will be set to the quantile $1-\alpha$ of the asymptotic distribution of Proposition \ref{prop:asymt} or \ref{prop:dw} or Corollary \ref{coro:dw}, depending on the context.

\section{Theory and Methodology} \label{Sec:Theory}
\subsection{Asymptotic Tests for Dimension}
Assume now that $\bo x$ satisfies Model \eqref{eq:model}. Then, let $q$ denote the true value of the signal dimension (i.e., $H_{0q}$ is true) and consider the limiting distribution of $t_q$. To establish the asymptotic results, we need to introduce a few technical conditions:

\begin{condition}\label{c:3}
The random fields $z_1,\cdots,z_p$ are independent centered stationary Gaussian random fields.
\end{condition}

The zero mean assumption is replaced by a constant unknown mean assumption in Section \ref{subsec:general}.
For $a=1,\cdots, p$, we let $z_a$ have stationary covariance function $K_a:\R^d\to\R$ with $\mathop{\mathrm{Cov}}\nolimits(z_a(\bo s), z_a(\bo s+\bo h))=K_a(\bo h).$

\begin{condition}\label{c:4}
A fixed $\delta>0$ exists such that, for all $n\in\mathbb{N}$ and, for all $i\neq j, i,j=1,\cdots,n, ||\bo s_i-\bo s_j||\ge\delta.$
\end{condition}

\begin{condition}\label{c:5}
Fixed $\beta>0$ and $\alpha>0$ exist such that, for all $\bo y\in\R^d$ and, for all $a=1,\cdots,p$
\[
|K_a(\bo y)|\le\frac{\beta}{1+||\bo y||^{d+\alpha}}.\]
\end{condition}

\begin{condition}\label{c:6}
Assuming Condition \ref{c:5} holds, then for the same $\beta>0$ and $\alpha>0$, we have, for $\ell=1,\cdots,k$
\[
|f_{\ell}(\bo y)|\le\frac{\beta}{1+||\bo y||^{d+\alpha}}.\]
\end{condition}

\begin{condition}\label{c:7}
For $\ell=1,\cdots,k$, we have
\[
\liminf_{n\to\infty}F_{n,f_{\ell}}>0.\]
\end{condition}

\begin{condition}\label{c:8}
For all $\ell, \ell^{'}=1,\cdots,k, \ell\neq \ell^{'}, f_{\ell}(\bo y)f_{\ell^{'}}(\bo y)=0$ for all $\bo y\in\R^d$.
\end{condition}

Condition \ref{c:4} implies that we are dealing with the increasing-domain asymptotic framework. For examples, see \cite{cressie2015statistics} for an introduction and \cite{bevilacqua2012estimating} for recent developments.
Condition \ref{c:5} holds for all the standard covariance functions in spatial statistics (\cite{cressie2015statistics}).
A typical example of a function $f\in\{f_1,\cdots,f_k\}$ for which Conditions \ref{c:2} and \ref{c:6} are satisfied is the ``ring" kernel:
\begin{equation}\label{eq:ring_kernel}
R(r_1, r_2)(\bo s)=I(r_1<||\bo s||\le r_2),
\end{equation}
with $0<r_1<r_2<\infty$. Condition \ref{c:7} is mild and simply requires that for $\ell=1,\cdots,k$, the number of pairs of observation locations $\bo s_i, \bo s_j, i, j=1,\cdots,n$, for which $f_{\ell}(\bo s_i-\bo s_j)$ is non-zero is not negligible when compared with $n$. Condition \ref{c:8} means that the supports of the kernels are disjoint. This enables us to have a simple and elegant asymptotic distribution of the test statistic.

Our first main result is on the asymptotic null distribution of our test statistic $t_q$.

\begin{proposition}\label{prop:asymt}
Assume that Conditions \ref{c:1}-\ref{c:8} hold. Then, as $n\to\infty$,
\[
t_q\xrightarrow{d}\chi^2_{k(p-q)(p-q+1)/2}.\]
\end{proposition}

In the next proposition, we show that when considering the same normalization as that considered by \cite{bachoc2018spatial} for the local covariance matrices, and when removing the assumption of disjoint kernel supports, we still obtain an asymptotic distribution of the test statistic as the distribution of the squared Euclidean norm of a Gaussian vector. In this proposition, we consider a metric $d_w$ generating the topology of weak convergence on the set of Borel probability measures on Euclidean spaces (e.g., \cite{dudley2018real}, p. 393).

\begin{proposition}\label{prop:dw}
Assume that Conditions \ref{c:1}-\ref{c:6} hold. Let the test statistic $\tilde{t}_q$ be defined as $t_q$, with $\widehat{\bo M}(f)$ replaced by
\[
\widetilde{\bo M}(f)=\frac{1}{n}\sum_{i=1}^n\sum_{j=1}^nf(\bo s_i-\bo s_j)\bo x(\bo s_i)\bo x(\bo s_j)^T\]
for $f\in\{f_1,\cdots,f_k\}$. Let $\mathcal{L}_{\tilde{t}_q, n}$ be the distribution of the test statistic $\tilde{t}_q$, and let $\mathcal{L}_{\bo V,n}$ be the distribution of
$
\sum_{\ell=1}^k\sum_{a,b=1}^{p-q}\bo V^2_{\ell,a,b},
$
where $(\bo V_{\ell,a,b})_{\ell=1,\cdots,k,a,b=1,\cdots,p-q}$ is a Gaussian vector with mean vector $\bm{0}$ and with covariance matrix defined by
\[
\mathop{\mathrm{Cov}}\nolimits(\bo V_{\ell,a,b}, \bo V_{\ell', a', b'})=\frac{1}{2}F_{n, f_{\ell}, f_{\ell'}}(I(a=a')I(b=b')+I(a=b')I(b=a'))\]
with
\[
F_{n, f_{\ell}, f_{\ell'}}=\frac{1}{n}\sum_{i,j=1}^nf_{\ell}(\bo s_i-\bo s_j)f_{\ell'}(\bo s_i-\bo s_j),\]
for $\ell, \ell'=1,\cdots,k$ and $a, b, a', b'=1,\cdots,p-q$. Then, as $n\to\infty$,
\[
d_w(\mathcal{L}_{\tilde{t}_q, n}, \mathcal{L}_{\bo V,n})\to 0.\]
\end{proposition}

In the following corollary, we show that if the supports of the kernels are disjoint, the test statistic converges to a weighted chi-squared distribution. See e.g. \cite{Bodenham-Adams:2016} for a presentation of the approximation procedures for this distribution.

\begin{corollary}\label{coro:dw}
Consider the setting of Proposition \ref{prop:dw} and assume additionally that Condition \ref{c:8} holds. Then, the limiting distribution $\mathcal{L}_{\bo V,n}$ in Proposition \ref{prop:dw} is equal to the distribution of
\[
\sum_{\ell=1}^kF_{n,f_{\ell}}\mathcal{X}_{\ell}^2\]
where $\mathcal{X}_1^2,\cdots,\mathcal{X}_k^2$ are independent and are chi-squared distributed with $(p-q)(p-q+1)/2$ degrees of freedom.
\end{corollary}

\subsection{Regular Domain as a Special Example}
When the data are observed in a regular-grid domain, i.e., $\mathcal{S}\subseteq\mathbb{Z}^d$, the kernel functions can be based on the natural notion of a spatial neighborhood on the grid, which simplifies our technique.

A location $\bo s_0=(s_1, \cdots, s_d)\in\mathbb{Z}^d$ has $2d$ one-way lag-$h$ neighbors, $(s_1\pm h, \cdots, s_d), (s_1, s_2\pm h, \cdots, s_d), \cdots, (s_1, \cdots, s_{d-1}, s_d\pm h)$. For example, if $d=2$ and $h=1$, the 4 one-way lag-$1$ neighbors of $\bo s_0$ are ``left" $(s_1-1, s_2)$, ``right" $(s_1+1, s_2)$, ``up" $(s_1, s_2+1)$ and ``down" $(s_1, s_2-1)$. Therefore, we could define the one-way lag-$1$ population and sample local covariance matrices as
\begin{small}
\begin{equation}
\bo M =
\frac{1}{\sqrt{n\sum_{i=1}^n|\mathcal{N}_{\bo s_i}|}}
\sum_{i=1}^n \sum_{\bo s_j\in \mathcal{N}_{\bo s_i}}
\mathbb{E}\left(\bo x(\bo s_i)\bo x(\bo s_j)^T\right)
\quad \text{and} \quad
\widehat{\bo M} =
\frac{1}{\sqrt{n\sum_{i=1}^n|\mathcal{N}_{\bo s_i}|}}
\sum_{i=1}^n\sum_{\bo s_j\in \mathcal{N}_{\bo s_i}}
\bo x(\bo s_i)  \bo x(\bo s_j)^T
\end{equation}
\end{small}
where, for $\bo x \in \mathbb{Z}^d$,
\[
\mathcal{N}_{\bo x}=\{\bo s \in \{ \bo s_1 , \ldots , \bo s_n \} ; |\bo x - \bo s | = 1 \},
\]
with $|\boldsymbol{u}| = |u_1| + \dots + |u_d|$ for $\boldsymbol{u} = (u_1,\ldots,u_d) \in \mathbb{R}^d$.
The matrices $\bo M$ and $\widehat{\bo M}$ are of the form $\bo M(f)$ and $\widehat{\bo M}(f)$ in \eqref{eq:plc} and \eqref{eq:slc} for $f(\bo s) = I( ||\bo s|| = 1 )$, $\bo s \in \mathbb{R}^d$.

Similarly, if $d=2$ a location $\bo s_0 = (s_1,s_2) \in \mathbb{Z}^2$ has $4$ two-way lag-$1$ neighbors that are of the form $(s_1 \pm 1 , s_2 \pm 1)$.
In general, for $m,h \in \mathbb{N}$, $1 \leq m \leq d$, the $m$-way lag-$h$ population and sample local covariance matrices can be defined as:
\begin{small}
\begin{equation}
\bo M=\frac{1}{\sqrt{n\sum_{i=1}^n|\mathcal{N}^m_{h, \bo s_i}|}}\sum_{i=1}^n\sum_{\bo s_j\in \mathcal{N}^m_{h, \bo s_i}}
\mathbb{E}\left(\bo x(\bo s_i)\bo x(\bo s_j)^T\right)
\quad \text{and} \quad
\widehat{\bo M}=\frac{1}{\sqrt{n\sum_{i=1}^n|\mathcal{N}^m_{h, \bo s_i}|}}\sum_{i=1}^n\sum_{\bo s_j\in \mathcal{N}^m_{h, \bo s_i}} \bo x(\bo s_i)\bo x(\bo s_j)^T,
\end{equation}
\end{small}
where, for $\bo x \in \mathbb{Z}^d$,
\[
\mathcal{N}^m_{h, \bo x}
=
\{
\bo s \in \{ \bo s_1 , \ldots , \bo s_n \};
 \bo s=\psi_J(\bo x, \zeta_J(\bo x) + h \boldsymbol{v}),
~ \text{for some} ~
  J \in \mathcal{A}_m, \boldsymbol{v} \in \{ - 1 , 1 \}^m
\},
\]
with $\mathcal{A}_m=\{J=(i_1,\cdots,i_m)\in\mathbb{N}^m; 1 \leq i_1<\cdots<i_m\le d\}$,
 that is
$|\mathcal{A}_m|=\binom{d}{m}$
and for $J =(i_1,\cdots,i_m) \in\mathcal{A}_m, \bo y=(y_1, \cdots, y_m)\in\mathbb{Z}^m, \zeta_J(\bo x)=(x_{i_1},\cdots,x_{i_m}), \psi_J(\bo x, \bo y)=(x_1, \cdots, x_{i_1-1}, y_1, x_{i_1+1}, \cdots, x_{i_m-1}, y_m, x_{i_m+1}, \cdots, x_d)$.

In general, in Equations~\eqref{eq:plc} and \eqref{eq:slc} the $m$-way lag-$h$ population and sample local covariance matrices $\bo M$ and $\widehat{\bo M}$ can also be written in the form  $\bo M(f)$ and $\widehat{\bo M}(f)$ , with $f(\bo s) = I( \bo s \in \{ -h,0,h \}^d, |\bo s| = hm )$, $\bo s \in \mathbb{R}^d$.

Consequently, for the similarly defined test statistic $t_q$, the same limiting conclusions can be derived. By exploiting the neighborhood structure in the regular domain case, we can also shorten the computation time for other techniques, such as the proposed asymptotic test or spatial bootstrap, with the greatest time improvement being achieved for the latter (see Sections \ref{subsec:bootstrap} and \ref{subsec:sim2}).

\subsection{Estimation of the Number of Signal Components}

In this section, we investigate an estimator of the signal number $q$ based on the asymptotic tests. Now we wish to test the null hypothesis, for $r\in\{0,\cdots,p-1\}$,
\begin{center}
$H_{0r}:$ There are exactly $p-r$ white noise processes in $\z$.
\end{center}
This hypothesis states that the signal dimension is $r$. Similar to Section \ref{Sec:Setup}, for $r = 0,\ldots,p-1$,  we can partition, for $\ell = 1 , \ldots , k$,
\[
\widehat \bo D_{\ell} =
\left(
\begin{array}{cc}
\widehat \bo D_{\ell,rr} & \widehat \bo D_{\ell,r-r}\\
\widehat \bo D_{\ell,-rr} & \widehat \bo D_{\ell,-r-r}\\
\end{array}
\right)
\]
where the block $\widehat{\bo D}_{\ell,rr}$ has size $r \times r$ and the block $\widehat \bo D_{\ell,-r-r}$ has size $(p - r) \times (p-r)$. Then, consider then the test statistic:
\[
t_r = \frac{n}{2} \sum_{\ell=1}^k ||\widehat \bo D_{\ell,-r-r}||^2.
\]

We now use the test statistic $t_r, r=0,1,\cdots, p-1$ for the estimation problem and derive a number of useful limiting properties in the following proposition.

\begin{proposition} \label{prop:power}
Assume the same conditions as in Proposition \ref{prop:asymt}. Then, \begin{itemize}
\item If $r \geq q$, then $t_r$ is bounded in probability.
\item If $r < q$, then there exists a fixed $b >0$ such that $ t_r / n \geq b + o_p(1)$.
\end{itemize}
\end{proposition}

A consistent estimate $\hat{q}$ of the unknown signal dimension $q\le p-1$ can then be based on the test statistic $t_r$ as follows.

\begin{proposition} \label{prop:estimating:q}
Assume the same conditions as in Proposition \ref{prop:asymt}.
Let $(c_n)_{n \in \mathbb{N}}$ be a sequence of positive numbers such that $c_n \to \infty$ and $c_n = o(n)$ as $n \to \infty$. Let
\[
\hat{q} = \min \left\{ r \in \{1,\dots,p-1\} \,|\,t_r \leq c_n \right\},
\]
 with the convention $\min \emptyset = p$. Then, $\hat{q} \to q$ in probability as $n \to \infty$.
\end{proposition}

Specifying the sequence $c_n$ is not obvious in practice.
Nevertheless, an estimator $\hat{q}$ can also be found by applying a suitable strategy to perform successive tests.
Later, in the simulations we will always test for simplicity at the same significance level and apply a divide-and-conquer strategy for the testing.

\subsection{General Mean}\label{subsec:general}

The previous results were derived under the assumption that $\mathbb{E}(\z(\bo s))=\bo 0$. In the next proposition, we show that the conclusions of Propositions \ref{prop:asymt}, \ref{prop:dw}, \ref{prop:power}, and \ref{prop:estimating:q} and of Corollary \ref{coro:dw} are unchanged when $\z$ has a non-zero unknown constant mean function and when the observations are empirically centered for the computation of the local covariance matrices.

\begin{proposition}\label{prop:nonzerom}
Assume that for $a=1,\cdots,p, z_a$ has constant mean function $\mu_a\in\R$. Let, for $f\in\{f_1,\cdots,f_k\},$
\begin{equation}\label{eq:Mbar}
\widebar{\bo M}(f)=\frac{1}{n\sqrt{F_{n,f}}}\sum_{i=1}^n\sum_{j=1}^nf(\bo s_i-\bo s_j)(\bo x(\bo s_i)-\bar{\bo x})(\bo x(\bo s_j)-\bar{\bo x})^T \quad \text{and} \quad \widebar{\bo M}(f_0)=\frac{1}{n}\sum_{i=1}^n(\bo x(\bo s_i)-\bar{\bo x})(\bo x(\bo s_i)-\bar{\bo x})^T,
\end{equation}
with $\bar{\bo x}=(1/n)\sum_{i=1}^n\bo x(\bo s_i)$.\\
Then, the conclusions of Propositions \ref{prop:asymt}, \ref{prop:dw}, \ref{prop:power}, and \ref{prop:estimating:q} as well as that of Corollary \ref{coro:dw} still hold under the same assumptions, except that $\widehat{\bo M}(f)$ is everywhere replaced by $\widebar{\bo M}(f)$.
\end{proposition}

For the remainder of the paper, we assume that the mean is unknown.

\subsection{Bootstrap Tests for Dimension}\label{subsec:bootstrap}

The above derived noise dimension test based on the large sample behavior of the introduced test statistic is efficient to compute, but a large sample size may be needed for the finite sample level to match the asymptotic one. As an alternative for smaller sample sizes, we can formulate noise dimension tests based on the bootstrap. These latter tests also have the benefit of not requiring $\z$ to be Gaussian.

In its original form, the bootstrap is a non-parametric tool for estimating the distribution of an estimator or test statistic by re-sampling from the empirical cumulative distribution function (ecdf) of the sample at hand. It has had good performance in many statistical problems by theoretical analysis as well as simulation studies and applications to real data. See \cite{chernick2011bootstrap} or \cite{lahiri} for a more detailed discussion.

Again, we assume that the observed random field is following the SBSS model given by Equation~\eqref{eq:model} and want to test $H_{0r}$ given an SBSS solution of Equation~\eqref{eq:hg} for a certain kernel setting and the corresponding test statistic seen in Equation~\eqref{eq:t}. In the following, we formulate a method for re-sampling from the distribution of Model \eqref{eq:model} by respecting the null hypothesis $H_{0r}$. In line with the ideas presented by \cite{bootstrap_ts} this is achieved by leaving the hypothetical signal part of the estimated latent field $\hat{\z}(\bo s) = \hat{ \G} \bo x(\bo s) $ untouched and manipulating only the hypothetical noise parts $(\hat{\z}(\bo s))_i$ for $i = r+1,\dots,p$ and all $\bo s \in \{\bo s_1 , \ldots , \bo s_n  \}$ in one of the following ways.

\textbf{Parametric:} Here, it is assumed that each noise part is independent and identically distributed (iid) Gaussian, as is usual for white noise processes. This leads to bootstrap samples $(\z^* (\bo s) )_{i} \sim N(0,1)$ for $i = r+1,\dots,p$ and corresponding to each $\bo s \in \{\bo s_1 , \ldots , \bo s_n  \}$.

\textbf{Permute:} Here, we assume that each noise component is still iid but that it does not necessarily follow a Gaussian distribution. Therefore, bootstrap samples are drawn from the ecdf of the joint noise components: $(\z^* (\bo s) )_{i} \sim \text{ecdf}((\hat{\z}(\bo s_1)^\top)_{\hat{w}}, \dots, (\hat{\z}(\bo s_n)^\top)_{\hat{w}})$, with $i=r+1,\ldots,p$, $\bo s \in \{ \bo s_1 , \ldots , \bo s_n \}$ and where $\hat{w}$ denotes the noise components ($r+1$ to $p$) of $\hat{\z}$.

\begin{algorithm}
	Set the number of resamples $B$, the observed sample $\bo X = (\bo x(\bo s_1), \ldots, \bo x(\bo s_n))^\top$, the flag $spatial \_ resampling$ and optionally the block size $m$;\;\\
	Compute the SBSS solution and get $\hat{\bo \Gamma}$ and $\hat{\bo Z} = (\hat{ \G} \bo X^\top)^\top $ and compute test statistic $t=t_r(\bo X)$;\;\\
	\For{$k \in \{ 1, \ldots , B \}$}{
	 Replace the last $p-r$ columns of $\hat{\bo Z}$ by either a parametric or bootstrap sample to get $\bo Z^{*k}$;\;\\
	\If{spatial\_resampling = TRUE}{Replace $\bo Z^{*k}$ by a full spatial bootstrap sample. See text for details.;}
	Compute $\bo X^{*k} \leftarrow \hat{\bo \Gamma} \bo Z^{*k} $ and $t^k \leftarrow t_r(\bo X^{*k})$;\;
}
	Return the $p$-value: $[\#(t^k \geq t) + 1]/(B + 1)$;\;
	\caption{Testing $H_{0r} : q=r$}
	\label{alg::boot}
\end{algorithm}

After replacing the hypothetical noise part by a bootstrap sample in one of the former ways, the goal of sampling from Model \eqref{eq:model} under $H_{0r}$ is achieved. However, so far the uncertainty of estimating the signal has not been considered in the bootstrap test. Therefore, an optional second step in the whole re-sampling procedure is devoted to drawing a spatial bootstrap sample from the already manipulated sample as follows. We suggest the application of spatial bootstrapping as discussed in \cite{lahiri}, and in the following we summarize the main ideas. Let us recall that the set of sampling sites $\mathcal{C} = \{ \bo s_1 , \ldots , \bo s_n \}$ lies inside the $d$-dimensional spatial domain $\mathcal{S}$, which can be viewed as the ``sample region'' and hence $\mathcal{C} \subseteq \mathcal{S} \subseteq \mathbb{R} ^ d$. $\mathcal{S}$ is divided into non-overlapping blocks of size $m^d$ that lie partially in $\mathcal{S}$, formally $\mathcal{B} = \{ b_{\bo i} =  (\bo i + (0, 1] ^ d)m \cap \mathcal{S} : (\bo i + (0, 1] ^ d)m \cap \mathcal{S} \neq \emptyset, \bo i \in \mathbb{Z}^d \}$, and overlapping blocks that lie fully in $\mathcal{S}$, written as $\mathcal{B}_{bs} = \{ b_{\bo j} = \bo j + (0, 1] ^ d m : \bo j + (0, 1] ^ d m \subseteq \mathcal{S}, \bo j \in \mathbb{Z}^d \}$. The bootstrapped spatial domain $\mathcal{S}^*$ is formed by replacing each block $b_{\bo i} \in \mathcal{B}$ with a randomly with replacement sampled block $b_{\bo j} \in \mathcal{B}_{bs}$ that is trimmed to the shape of $b_{\bo i}$ by $b_{\bo j} \cap (b_{\bo i} - \bo i m + \bo j)$. Hence, the trimmed version of $b_{\bo j}$ remains at the original location of $b_{\bo j}$, while the shape changes to that of $b_{\bo i}$, taking care of the boundary blocks that do not fully lie within $\mathcal{S}$. Finally, the bootstrapped version of the random field writes as $\bo z^* = \{\bo z(\bo s) : \bo s \in \mathcal{S}^* \cap \mathcal{C} \}$. Note that in each spatial bootstrap iteration, the shape of $\mathcal{S}^*$ and therefore the bootstrapped sampling sites differ. This in turn makes the computation of the local covariance matrices a demanding task, as it relies on the distances between all sampling sites, which need to be newly computed in each iteration. For regular data, this can be avoided by using a slightly different bootstrap regime as follows.

\cite{lahiri_regular} have suggested a slightly different approach for sampling sites located on a regular grid, meaning that the sampling sites satisfy $\{ \bo s_1, \dots, \bo s_n \} \subseteq \mathcal{S} \cap \mathbb{Z}^d$. Again, the domain $\mathcal{S}$ is divided into blocks of size $m ^ d$ that are either non-overlapping or overlapping but lie completely inside $\mathcal{S}$, leading to $\mathcal{B} = \{ (\bo i + (0, 1] ^ d)m : (\bo i + (0, 1] ^ d)m \subseteq \mathcal{S}, \bo i \in \mathbb{Z}^d \}$ and $\mathcal{B}_{bs}$, as defined above. The key difference is that the bootstrap sample is drawn at the level of the random field values, whereas the former bootstrap version operates at the level of the spatial domain. Specifically, for each block $b_{\bo i} \in \mathcal{B}$ the values $\{\z(\bo s) : \bo s \in b_{\bo i} \cap \mathbb{Z}^d\}$ are replaced by $\{\z(\bo s) : \bo s \in b_{\bo j} \cap \mathbb{Z}^d \}$ for a randomly with replacement chosen block $b_{\bo j} \in \mathcal{B}_{bs}$. This procedure keeps the bootstrapped spatial domain and sampling sites equal in all iterations, namely the unison of all blocks from $\mathcal{B}$. This in turn simplifies the computation of local covariance matrices, as only the random field values change. We will compare the computation times of the former two approaches in the simulation study presented in Section~\ref{subsec:sim2}.

Algorithm~\ref{alg::boot} summarizes the formerly discussed bootstrap strategy to test for one specific value of signal dimension $r$. To estimate the signal dimension, a sequence of tests for different signal dimensions $r$ at a given significance level $\alpha$ are carried out. A number of different test sequences are possible, but we rely on a divide-and-conquer strategy outlined in Algorithm~\ref{alg::div_con}. Here, the $test \_ function$ could be either one of the bootstrap test variants seen in Algorithm~\ref{alg::boot} or the asymptotic test outlined above.

\begin{algorithm}
	Set $lower$, $upper$ and $\alpha$;\;\\
	$middle = \lfloor (upper-lower)/2 \rfloor$;\;\\
	\While{$(middle != lower) ~ \& \& ~ (middle != upper)$}{
	 $p = test \_ function(r=middle)$;\;\\
	\If{$p < alpha$}{$lower = middle$;}
	\Else{$upper = middle$;}
	$middle = \lfloor (upper-lower)/2 \rfloor$;\;	}
	Return $\hat{q} = middle + 1$;\;
	\caption{Divide and Conquer}
	\label{alg::div_con}
\end{algorithm}

\section{Simulation}\label{Sec:Sim}

To validate the performance of the methods introduced above, we carried out three extensive simulation studies in R 3.6.1 (\cite{r_language}) with the help of the packages SpatialBSS (\cite{spatialbss_package}), JADE (\cite{jade_package}), sp (\cite{sp_package}), raster (\cite{raster_package}), gstat (\cite{gstat_package}) and RandomFields (\cite{randomfields_package}).

\subsection{Simulation Study 1: Hypothesis Testing}\label{sec::sim_hypothesis}
In this part of the simulation, we explored the performance of hypothesis testing. For all the following simulations, we considered the SBSS model, as shown in Equation~\eqref{eq:model}, where without loss of generality we set $\mu_a = 0$ for $a=1,\dots,p$ and assume the mean to be unknown. For the latent signal part we used two different three-variate random field model settings. Therefore, the true dimension is always $q=3$. All the random fields followed a Mat\'ern correlation structure, and the $a$-th random field $z_a$ thus had its covariance function value at $\boldsymbol{u}, \boldsymbol{v} \in \mathbb{R}^d$, given by:
\[
K_a(h; \nu, \phi) = \frac{1}{2 ^ {\nu - 1} \Gamma (\nu)} \left( \frac{h}{\phi} \right) ^ \nu  K_\nu \left( \frac{h}{\phi} \right),
~ ~
h = || \boldsymbol{u} - \boldsymbol{v} ||,
\]
where $\nu >0$ is the shape parameter, $\phi >0$ is the range parameter, and $K_\nu$ is the modified Bessel function of second kind with shape parameter $\nu$. The parameters used were $(\nu, \phi) \in \{(3,2), (2,1.5), (1,1) \}$ and $\{(3,2), (2,1.5), (0.6,0.6) \}$ for model setting 1 and 2 respectively, which are depicted in Figure~\ref{fig:coordinates}. Model setting 2 can be viewed as a low-dependence version of model setting 1. The noise part always consists of iid samples drawn from $N_{2}(\bo 0, \I_{2})$, leading to a total latent field dimension of $p=5$ for both model settings. As SBSS is affine equivariant (for details see \cite{bachoc2018spatial} and the appendix) we chose the mixing matrix to be the identity matrix, i.e., $\Om = \bf I_5$, without loss of generality.

We focused on the squared spatial domains $[0,n] \times [0,n]$ (also written in the following as $n \times n$) of different sizes $n \in \{30, 40, 50, 60\}$. For a given domain, we considered two different sample location patterns: uniform and skewed. For the uniform pattern, $n ^ 2$ pairs of $(x,y)$-coordinates were randomly drawn from a uniform distribution $U(0,1)$ and then multiplied by $n$, leading to a constant sampling location density over the entire domain. We followed the same approach for the skewed pattern, with the only difference being that the $x$ coordinate values were drawn from a beta distribution $\beta (2,5)$, resulting in a denser arrangement of samples in the left half of the domain.

For the local covariance matrices \eqref{eq:slc}, we used two different kernel function settings. Kernel setting 1 used only one ring kernel function \eqref{eq:ring_kernel} with parameters $(r_1,r_2) = (0,2)$, while kernel setting 2 used three ring kernel functions with parameters $(r_1,r_2) \in \{(0, 2), (2, 4), (4, 6)\}$. Figure~\ref{fig:coordinates} depicts a simulation example for each of the uniform and skewed coordinate patterns, where the circles represent the different ring kernel radii.

\begin{figure}
\centering
    \begin{minipage}[t]{0.32\linewidth}
        \centering
        \includegraphics[scale=0.5]{figures/cov_models.pdf}
    \end{minipage}
    \hfill
    \begin{minipage}[t]{0.32\linewidth}
        \centering
        \includegraphics[scale=0.5]{figures/unif.pdf}
    \end{minipage}
    \hfill
    \begin{minipage}[t]{0.32\linewidth}
        \centering
        \includegraphics[scale=0.5]{figures/skew.pdf}
    \end{minipage}
    \caption{
Left: Mat\'ern correlation functions for model setting 1, which consists of the signal random field $(z_1,z_2,z_{3,1})$ with parameters $(\nu, \phi) \in \{ (3,2), (2,1.5), (1,1) \}$ and model setting 2 formed by the signal random field $(z_1,z_2,z_{3,2})$ with parameters $(\nu, \phi) \in \{(3,2), (2,1.5), (0.6,0.6)\}$.
Middle and right: uniform (middle) and skewed (right) coordinate sample pattern for a spatial domain of size $30 \times 30$ with three circles of radii $(2,4,6)$ representing ring kernel functions.}
        \label{fig:coordinates}
\end{figure}

For each of the four simulation settings, we carried out 2000 repetitions, and in each repetition we tested three different null hypothesis ($H_{02}$, $H_{03}$, and $H_{04}$) with the following five test approaches: asymptotic test (Asym), noise bootstrapping with option parametric (Param), noise bootstrapping with option permute (Perm), full spatial bootstrapping with option parametric (Sp Param), and full spatial bootstrapping with option permute (Sp Perm). For all bootstrap approaches, we fixed the number of re-samples to be $B=200$, and for the full spatial bootstrap, the block size was equal to $m=10$.

Rejection rates based on a significance level of $\alpha = 0.05$ for all simulation settings are presented in Tables~\ref{tab:model_1} and \ref{tab:model_2}. Overall, all the test methods appeared to maintain the expected rejection rates, which were $1.00$ for $H_{02}$, $0.05$ for $H_{03}$, and $<0.05$ for $H_{04}$ based on $\alpha = 0.05$. Only for small samples sizes ($30 \times 30$) did the asymptotic test show a too small rejection rate for kernel setting 2 and the skewed sample location pattern. Thus, for practical applications, smaller numbers of kernel functions might be preferable for the asymptotic test. For bootstrapping, the full spatial variants and those relying only on manipulating the hypothetical noise part performed equally well. Considering the computation time, the latter bootstrap variant might be preferable, as explored in more detail in Section \ref{subsec:sim2}.

\begin{table}[h]
\caption{Rejection rates for model setting 1  at level $\alpha=0.05$.}
\label{tab:model_1}
\centering
\scalebox{0.8}{
\begin{tabular}{clcccccccccccc}
  \toprule
  &  &  \multicolumn{6}{c}{Uniform} & \multicolumn{6}{c}{Skew} \\
   \cmidrule(r){3-8} \cmidrule(l){9-14}
   &  & \multicolumn{3}{c}{Kernel Setting 1} & \multicolumn{3}{c}{Kernel Setting 2} & \multicolumn{3}{c}{Kernel Setting 1} & \multicolumn{3}{c}{Kernel Setting 2}\\
   \cmidrule(r){3-5} \cmidrule(l){6-8} \cmidrule(l){9-11} \cmidrule(l){12-14}
Domain & Method & $H_{02}$ & $H_{03}$ & $H_{04}$ & $H_{02}$ & $H_{03}$  & $H_{04}$ & $H_{02}$ & $H_{03}$ & $H_{04}$ & $H_{02}$ & $H_{03}$  & $H_{04}$ \\
  \midrule
\multirow{5}{*}{$30 \times 30$} & Asym & 1.000 & 0.041 & 0.006 & 1.000 & 0.042 & 0.007  & 1.000 & 0.042 & 0.004 & 1.000 & 0.029 & 0.003\\
   & Sp Param & 1.000 & 0.048 & 0.006 & 1.000 & 0.058 & 0.001 & 1.000 & 0.059 & 0.004 & 1.000 & 0.051 & 0.000\\
   & Sp Perm & 1.000 & 0.050 & 0.006 & 1.000 & 0.059 & 0.000 & 1.000 & 0.058 & 0.004 & 1.000 & 0.052 & 0.001\\
   & Param & 1.000 & 0.042 & 0.006 & 1.000 & 0.044 & 0.006 & 1.000 & 0.050 & 0.008 & 1.000 & 0.039 & 0.005\\
   & Perm & 1.000 & 0.045 & 0.008 & 1.000 & 0.051 & 0.006 & 1.000 & 0.049 & 0.008 & 1.000 & 0.035 & 0.005\\ \midrule
  \multirow{5}{*}{$40 \times 40$} & Asym & 1.000 & 0.055 & 0.004 & 1.000 & 0.048 & 0.005 & 1.000 & 0.045 & 0.002 & 1.000 & 0.040 & 0.005\\
   & Sp Param & 1.000 & 0.056 & 0.005 & 1.000 & 0.066 & 0.000 & 1.000 & 0.056 & 0.003 & 1.000 & 0.064 & 0.002\\
   & Sp Perm & 1.000 & 0.063 & 0.005 & 1.000 & 0.061 & 0.000 & 1.000 & 0.055 & 0.004 & 1.000 & 0.065 & 0.002 \\
   & Param & 1.000 & 0.052 & 0.007 & 1.000 & 0.055 & 0.003 & 1.000 & 0.050 & 0.007 & 1.000 & 0.048 & 0.005\\
   & Perm & 1.000 & 0.056 & 0.007 & 1.000 & 0.052 & 0.004 & 1.000 & 0.048 & 0.008 & 1.000 & 0.050 & 0.004\\ \midrule
  \multirow{5}{*}{$50 \times 50$} & Asym & 1.000 & 0.049 & 0.005 & 1.000 & 0.040 & 0.010 & 1.000 & 0.040 & 0.006 & 1.000 & 0.044 & 0.009\\
   & Sp Param & 1.000 & 0.052 & 0.004 & 1.000 & 0.053 & 0.002 & 1.000 & 0.047 & 0.006 & 1.000 & 0.064 & 0.002\\
   & Sp Perm & 1.000 & 0.050 & 0.005 & 1.000 & 0.053 & 0.002 & 1.000 & 0.045 & 0.005 & 1.000 & 0.061 & 0.002\\
   & Param & 1.000 & 0.052 & 0.007 & 1.000 & 0.049 & 0.007 & 1.000 & 0.042 & 0.007 & 1.000 & 0.054 & 0.007\\
   & Perm & 1.000 & 0.050 & 0.008 & 1.000 & 0.050 & 0.006  & 1.000 & 0.042 & 0.010 & 1.000 & 0.054 & 0.008\\ \midrule
  \multirow{5}{*}{$60 \times 60$} & Asym & 1.000 & 0.052 & 0.006 & 1.000 & 0.048 & 0.010 & 1.000 & 0.044 & 0.004 & 1.000 & 0.045 & 0.004\\
   & Sp Param & 1.000 & 0.056 & 0.006 & 1.000 & 0.058 & 0.003 & 1.000 & 0.048 & 0.005 & 1.000 & 0.060 & 0.000\\
   & Sp Perm & 1.000 & 0.055 & 0.007 & 1.000 & 0.057 & 0.002 & 1.000 & 0.052 & 0.004 & 1.000 & 0.058 & 0.000\\
   & Param & 1.000 & 0.049 & 0.009 & 1.000 & 0.054 & 0.006 & 1.000 & 0.043 & 0.006 & 1.000 & 0.048 & 0.004\\
   & Perm & 1.000 & 0.053 & 0.009 & 1.000 & 0.050 & 0.008 & 1.000 & 0.046 & 0.006 & 1.000 & 0.048 & 0.004\\
  \bottomrule
\end{tabular}
}
\end{table}

\begin{table}[h]
\caption{Rejection rates for model setting 2 at level $\alpha=0.05$.}
\centering
\label{tab:model_2}
\scalebox{0.8}{
\begin{tabular}{clcccccccccccc}
  \toprule
    &  &  \multicolumn{6}{c}{Uniform} & \multicolumn{6}{c}{Skew} \\
     \cmidrule(r){3-8} \cmidrule(l){9-14}
   &  & \multicolumn{3}{c}{Kernel Setting 1} & \multicolumn{3}{c}{Kernel Setting 2} & \multicolumn{3}{c}{Kernel Setting 1} & \multicolumn{3}{c}{Kernel Setting 2}\\
   \cmidrule(r){3-5} \cmidrule(l){6-8} \cmidrule(l){9-11} \cmidrule(l){12-14}
Domain & Method & $H_{02}$ & $H_{03}$ & $H_{04}$ & $H_{02}$ & $H_{03}$  & $H_{04}$ & $H_{02}$ & $H_{03}$ & $H_{04}$ & $H_{02}$ & $H_{03}$  & $H_{04}$ \\
  \midrule
\multirow{5}{*}{$30 \times 30$} & Asym & 1.000 & 0.051 & 0.005 & 1.000 & 0.052 & 0.004 & 1.000 & 0.048 & 0.006 & 1.000 & 0.033 & 0.003\\
   & Sp Param & 1.000 & 0.053 & 0.005 & 1.000 & 0.062 & 0.000 & 1.000 & 0.058 & 0.005 & 1.000 & 0.055 & 0.002\\
   & Sp Perm & 1.000 & 0.052 & 0.006 & 1.000 & 0.065 & 0.001 & 1.000 & 0.056 & 0.006 & 1.000 & 0.051 & 0.001\\
   & Param & 1.000 & 0.052 & 0.011 & 1.000 & 0.058 & 0.003 & 1.000 & 0.059 & 0.011 & 1.000 & 0.043 & 0.004\\
   & Perm & 1.000 & 0.048 & 0.011 & 1.000 & 0.060 & 0.002 & 1.000 & 0.061 & 0.012 & 1.000 & 0.044 & 0.003\\ \midrule
  \multirow{5}{*}{$40 \times 40$} & Asym & 1.000 & 0.060 & 0.004 & 1.000 & 0.052 & 0.005 & 1.000 & 0.050 & 0.004 & 1.000 & 0.038 & 0.007\\
   & Sp Param & 1.000 & 0.063 & 0.002 & 1.000 & 0.060 & 0.000 & 1.000 & 0.060 & 0.004 & 1.000 & 0.054 & 0.002\\
   & Sp Perm & 1.000 & 0.055 & 0.002 & 1.000 & 0.062 & 0.000 & 1.000 & 0.058 & 0.002 & 1.000 & 0.057 & 0.002\\
   & Param & 1.000 & 0.056 & 0.006 & 1.000 & 0.056 & 0.004 & 1.000 & 0.052 & 0.008 & 1.000 & 0.045 & 0.005\\
   & Perm & 1.000 & 0.058 & 0.005 & 1.000 & 0.053 & 0.004 & 1.000 & 0.054 & 0.006 & 1.000 & 0.045 & 0.005\\ \midrule
  \multirow{5}{*}{$50 \times 50$} & Asym & 1.000 & 0.045 & 0.004 & 1.000 & 0.047 & 0.004 & 1.000 & 0.044 & 0.005 & 1.000 & 0.044 & 0.004\\
   & Sp Param & 1.000 & 0.048 & 0.002 & 1.000 & 0.056 & 0.000 & 1.000 & 0.053 & 0.002 & 1.000 & 0.058 & 0.001\\
   & Sp Perm & 1.000 & 0.049 & 0.002 & 1.000 & 0.053 & 0.001 & 1.000 & 0.050 & 0.005 & 1.000 & 0.055 & 0.001\\
   & Param & 1.000 & 0.045 & 0.004 & 1.000 & 0.050 & 0.002 & 1.000 & 0.048 & 0.007 & 1.000 & 0.051 & 0.004\\
   & Perm & 1.000 & 0.044 & 0.007 & 1.000 & 0.048 & 0.003 & 1.000 & 0.046 & 0.009 & 1.000 & 0.052 & 0.004\\ \midrule
  \multirow{5}{*}{$60 \times 60$} & Asym & 1.000 & 0.048 & 0.004 & 1.000 & 0.059 & 0.008 & 1.000 & 0.047 & 0.004 & 1.000 & 0.042 & 0.006\\
   & Sp Param & 1.000 & 0.052 & 0.005 & 1.000 & 0.072 & 0.002 & 1.000 & 0.050 & 0.004 & 1.000 & 0.059 & 0.000\\
   & Sp Perm & 1.000 & 0.056 & 0.003 & 1.000 & 0.068 & 0.002 & 1.000 & 0.050 & 0.004 & 1.000 & 0.057 & 0.000\\
   & Param & 1.000 & 0.047 & 0.009 & 1.000 & 0.063 & 0.004 & 1.000 & 0.046 & 0.005 & 1.000 & 0.052 & 0.003\\
   & Perm & 1.000 & 0.048 & 0.010 & 1.000 & 0.063 & 0.006 & 1.000 & 0.048 & 0.005 & 1.000 & 0.050 & 0.005\\
   \bottomrule
\end{tabular}
}
\end{table}

\subsection{Simulation Study 2: Computation Time Comparison}\label{subsec:sim2}

In this simulation, we investigated the computation times for the various test methods. As an illustrative example, we again considered a five-variate latent random field with model setting 1 and bivariate Gaussian noise components. In addition, we kept the same spatial domain sizes, though the sampling sites were changed to be regular defined as $[0,n] \times [0,n] \cap \mathbb{Z}^2$. $H_{03}$ was tested using the five former mentioned test methods with the same number of bootstrap samples and block sizes. The key difference is that each test was carried out with code designed for irregular sample locations as well as code that takes into account simplifications made possible by the fact that the sample locations were regular (e.g., the simplified spatial bootstrap algorithm). Two ring kernel functions with parameters $(r_1, r_2) \in \{(0, 1), (1, \sqrt{2}) \}$ were considered for the irregular code, and kernels of the form $f(\bo s)=I(||\bo s||=h)$ with $h \in \{1, \sqrt{2} \}$ were considered for the regular code (one-way and two-way lag-$1$ local covariance matrices). This choice ensured that the same neighbors were selected for both versions of the code and thus that the qualitative results of the tests were equal up to random effects of the bootstrap sampling procedures.

\begin{figure}[th]
	\centering
  \includegraphics[scale=0.5]{figures/running_times_2_k.pdf}
  \caption{Median running times of the five different test methods for different domain sizes with regular sampling sites based on five simulation repetitions. Computations were carried out with code designed for regular and irregular sampling sites.}
  \label{fig::2k_running_times}
\end{figure}

Figure~\ref{fig::2k_running_times} shows the median computation time based on five simulation repetitions carried out on a Windows machine with an Intel i5 CPU. The computation times revealed that asymptotic tests are fastest, as the SBSS solution needs to be computed only once, whereas bootstrap algorithms compute the SBSS solution $B$ times.

Of greater interest is the overall difference in the computation time between regular and irregular code. This might be explained by the fact that the code for regular sampling sites does not rely on distances between sampling sites as the irregular code does. Specifically, the selection of neighbors for local covariance matrices can be implemented by shifting the coordinate system appropriately for the regular code, whereas in the irregular code this is based on looping over the distance matrix among all coordinates. This difference should also explain the different scaling of the computation time with increasing sample size, as looping through the distance matrix depends on the actual number of locations, while coordinate shifting does not.

Further, there was a larger computation time difference between the full spatial bootstrap and the one that manipulates only the hypothetical noise for the irregular code compared with the regular one. This might be the impact of the simplified spatial bootstrap variant for regular sampling sites. As explained above, for the irregular code the distance matrix has to be computed again for every new iteration because the spatial bootstrap changes sampling sites for each iteration, which is not the case for the regular code, for which the sampling sites remain equal for each bootstrap iteration.

Overall, this simulation strongly indicates that regular sampling sites should be computationally treated as such. In addition, considering the overall similar performances of the tests in the former simulation, the spatial bootstrapping step for the irregular data might be discarded, as it significantly increases the computation time.

\subsection{Simulation Study 3: Estimation of the Signal Dimension}\label{subsec:sim3}

The former simulations investigated only hypothesis tests for one specific value of the hypothetical signal dimension. In this section, we explore the use of hypothesis tests for signal dimension estimation. We considered the exact same simulation settings as in Section~\ref{sec::sim_hypothesis} but increased the dimension of the noise part to seven leading to a total latent random field dimension of $p=10$, while the true signal dimension remained $q=3$. Estimation of the signal dimension was based on the divide-and-conquer strategy described above. All hypothesis tests were carried out using the asymptotic test method and the parametric bootstrap without full spatial bootstrapping. This choice is justified by the similar performance in signal dimension testing of all bootstrap test variants and the fact that the full spatial bootstrap is computationally unfeasible for such a large simulation.

\begin{figure}[ht]
	\centering
  \includegraphics[scale=0.5]{figures/estimation_1.pdf}
  \caption{Frequencies of the estimated signal dimension for model setting 1.}
  \label{fig::est_1}
\end{figure}

\begin{figure}[ht]
	\centering
  \includegraphics[scale=0.5]{figures/estimation_2.pdf}
  \caption{Frequencies of the estimated signal dimension for model setting 2.}
  \label{fig::est_2}
\end{figure}

Figures~\ref{fig::est_1} and \ref{fig::est_2} depict the estimated dimensions for 2000 simulation repetitions for a significance level of $\alpha=0.05$.
Overall, the estimation was highly accurate, with the estimated dimension being equal to the true one in approximately $95 \%$ of the cases.
 Interestingly, the signal dimension was never underestimated, while it was overestimated in approximately 100 of the simulation iterations, reflecting the significance level $\alpha=0.05$. For all settings, the asymptotic test performed better than the bootstrap test. This was especially true for low sample sizes, which is a counterintuitive result. However, it may be due to the fact that, as the former simulations show, for low sample sizes the asymptotic test never met the theoretical rejection rate, which is simply the significance level when the null is actually true for small sample sizes (Tables~\ref{tab:model_1} and \ref{tab:model_2}). Therefore, the true null is more often accepted leading to a better performance when estimating the signal dimension.

\section{Data Example}\label{Sec:Ex}

In this section, we illustrate the application of our tests to a geochemical dataset. Specifically, we analyze samples of the amount of 31 elements in terrestrial moss, which were collected during the Kola project at 594 sites located in Norway, Finland, and Russia alongside the coast of the Barents Sea. The data is freely available in the R package \cite{statda_package} and is described in more detail by \cite{moss_book}. \cite{nordhausen2015blind} and \cite{bachoc2018spatial} have already considered the dataset in the context of SBSS where the goal of the former publication was to identify geochemical interpretative components. Indeed, six meaningful components were found based on an SBSS solution that used a single local covariance matrix with a ball kernel $B(\bo s) = I(||\bo s|| \leq 50 \text{km})$, chosen based on experts knowledge. The latter publication introduces a variety of kernel function options and expands SBSS for the use of more than one local covariance matrix. The authors also found that jointly diagonalizing more than one local covariance matrix yielded a more stable solution compared with using only one. For the moss dataset, an SBSS solution based on four local covariance matrices exclusively using ring kernel functions with parameters $(r_1,r_2) \in \{(0, 25), (25, 50), (50, 75), (75, 100)\}$ km was used, among others, where the first six latent field components showed high correlation with the one presented by \cite{nordhausen2015blind}.

This was the starting point of our analysis, in which we estimated the signal dimension using the five test methods defined in Section \ref{sec::sim_hypothesis} based on the two kernel function settings described above. Note that, theoretically, the ball kernel function is not compatible with Condition~\ref{c:2}. Hence, we used exclusively ring kernel functions, leading to the two considered kernel settings: $(r_1,r_2)= (0, 25)$ (Setting 1) and $(r_1,r_2) \in \{(0, 25), (25, 50), (50, 75), (75, 100)\}$ (Setting 2) in km. Using the same methods as \cite{nordhausen2015blind} and \cite{bachoc2018spatial}, we pre-processed the data by performing an isometric log-ratio (ilr) transformation using pivot coordinates to respect the compositional nature of the geochemical data, which reduced the dimension of the data to $p=30$. The details of the ilr transformation and on compositional data in general can be found in \cite{comp_book}. For the bootstrap methods $B=200$ re-samples were drawn, and for the spatial bootstrap the domain was overlaid by a square grid with a side-length of 30 km. The block size was chosen to be $m=60$ km, meaning that every 30 km, a block of size 60 km was placed, forming the set of blocks from which the samples were drawn. All lengths were referred to the UTM zone 35N.

Table \ref{tab:kola} summarizes the estimated signal dimensions for the various test methods and kernel function settings based on a significance level of $\alpha = 0.05$. For both kernel settings, the dimension of the actual signal was approximately half of the original data dimension. This is a significant reduction in dimension, as further analysis (such as spatial prediction) needs to consider drastically fewer signal fields. Still, the estimated number of signals was higher than the number of meaningful geochemical components found by \cite{nordhausen2015blind}. The estimated number of signal dimensions for kernel setting 1 was always $\hat{q} = 14$, meaning it appears to be more stable (less sensitive to the choice of the test version) than estimations based on a greater number of ring kernel functions (setting 2), where the number ranged between $\hat{q} = 15$ and $18$. Thus, the results indicate that, overall, signal dimension estimations based on SBSS solutions with less kernel functions are more stable. This was already hinted in the former simulations and contrasts with the guidelines presented by \cite{bachoc2018spatial} where it was shown that SBSS solutions based on jointly diagonalizing more than one kernel function yielded more stable solutions.

\begin{table}
\caption{Estimated signal dimension for the moss dataset.}
\centering
\begin{tabular}{cccccc}
  \toprule
Kernel setting & Asym & Perm & Param & Sp Perm & Sp Param  \\ \midrule
Setting 1 & 14 & 14 & 14 & 14 & 14 \\
Setting 2 & 17 & 17 & 18 & 15 & 16 \\
  \bottomrule
\end{tabular}
\label{tab:kola}
\end{table}

Figures~\ref{fig::kola_1} and \ref{fig::kola_2} show the entries of the latent fields, and Figures~\ref{fig::kola_vars_1} and \ref{fig::kola_vars_2} show the corresponding sample variograms at the change point between the estimated signal and noise for both kernel function settings. Based on visual inspection of Figure~\ref{fig::kola_1}, IC.14-IC.16 might show some weak spatial dependence. However, the sample variograms in Figure~\ref{fig::kola_vars_1} suggest that IC.14 and IC.15 carry a very weak signal, and IC.16 already shows very similar behavior to the last part of the noise, IC.30. Each panel in Figure~\ref{fig::kola_2} could indicate very similar behavior and also illustrates the fact that different signal dimensions are estimated by the different methods. Interestingly, the sample variograms in Figure~\ref{fig::kola_vars_2} show that among IC.15-IC.20, IC.19 had the highest spatial dependence, though all seem to be highly similar to the last components.

\begin{figure}[h]
	\centering
  \includegraphics[scale=0.6]{figures/kola_1.pdf}
  \caption{Latent field components of the change-point between estimated signal and noise components for kernel function setting 1. Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.}
  \label{fig::kola_1}
\end{figure}

\begin{figure}[h]
	\centering
  \includegraphics[scale=0.6]{figures/variograms_1.pdf}
  \caption{Sample variograms for the latent field components that are first, around the change-point between the estimated signal and noise components and last for kernel function setting 1.}
  \label{fig::kola_vars_1}
\end{figure}

\begin{figure}[h]
	\centering
  \includegraphics[scale=0.6]{figures/kola_2.pdf}
  \caption{Latent field components of the change-point between the estimated signal and noise components for kernel function setting 2. Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.}
  \label{fig::kola_2}
\end{figure}

\begin{figure}[h]
	\centering
  \includegraphics[scale=0.6]{figures/variograms_2.pdf}
  \caption{Sample variograms for the latent field components that are first, around the change-point between the estimated signal and noise components and last for kernel function setting 2.}
  \label{fig::kola_vars_2}
\end{figure}

\section{Concluding Remarks}\label{Sec:Conclusion}
In this paper, we propose and study testing and estimation methods for the number of latent signal components in the SBSS model. The asymptotic null distributions of the test statistic are given under various conditions without assuming the domain is necessarily regular. A consistent estimator of the dimension based on the sequential tests is also introduced. For small sample cases, different bootstrap strategies are suggested. Besides the theoretical results, the three simulation studies presented in Section \ref{Sec:Sim} demonstrate that our asymptotic tests are comparable to the bootstrap ones in terms of hypothesis testing and estimation. In terms of computation time, our asymptotic method is much faster than the bootstrap ones. When a regular domain structure is used, the computation time can be significantly decreased.

Our proposed dimension tests in the SBSS context might be very useful for further analysis of the latent fields,  including for various forms of spatial prediction. Indeed, the components of the latent field are uncorrelated, and thus predictions can be carried out on each latent field independently, leading to a reduction from building a single multivariate model to building several univariate ones. This procedure was already investigated and found to be useful by \cite{kriging_paper}. As an additional step, one of our proposed dimension tests can be carried out before the spatial prediction, leading to a reduction of the latent field dimension, which results in the need for even fewer univariate models to be built.

Finally, in further research we plan to develop a ladle estimator (\cite{LuoLi:2016}) for this setting, as well as to develop similar approaches for spatiotemporal data. We will also study the fixed-domain asymptotic properties (\cite{cressie2015statistics}) of SBSS.

\section{Acknowledgments}
The work of CM, KN, and MY was supported by the Austrian Science Fund P31881-N32.