text:
\section*{Acknowledgment}
The authors would like to thank dr. Peyman Mohajerin Esfahani of the Delft Center for Systems and Control at the Delft University of Technology for his insightful comments.

\section{Theoretical analysis of \DBay}
\label{section:analysis}

\structure{General introduction}
This section will analyse the convergence of \DBay to the global optimum of a continuous DCOP in two parts.
Firstly, in \cref{subsection:convergenceofBO}, the convergence to the global optimum of the utility functions within the sampling phase is proven.
It will be shown that if the Lipschitz constant of the utility functions is known, the convergence to the global optimum can be guaranteed through the appropriate selection of a kernel and an acquisition function.
In this work, all utility functions are assumed to be Lipschitz continuous with known Lipschitz constant.
A utility function \mathcmd{\utilityFunction(\cdot)} is Lipschitz continuous with Lipschitz constant \lipschitzConstantIndexed if
\begin{align}
    \lipschitzConstantIndexedLong \qquad \forall \mathcmd{\functionInputBase_\mathcmd{i}}, \mathcmd{\functionInputBase_\mathcmd{j}} \in \domainOperator(\mathcmd{f})
    \label{eq:LipschitzConstant}
\end{align}
where $\domainOperator(\mathcmd{f})$ denotes the domain of the utility function.

Secondly, in \cref{subsection:globalconvergenceofDBAY}, the convergence of \DBay to the global optimum of the objective function based on the global optima of the utility functions is proven.
This analysis focusses on the assignment propagation phase.
The two parts of the analysis are combined to prove convergence of \DBay to the global optimum of continuous DCOPs with utility functions with known Lipschitz constants.

\input{sections/Analysis/section_convergence_global_utility_function.tex}

\input{sections/Analysis/section_convergence_global_objective_function.tex}

\subsection{Summary}
\structure{conclusion for a f with L we can select the kernel with parameter to guarantee global convergence!}
In \cref{subsection:convergenceofBO} it was shown that, based on the Lipschitz constant of a utility function (or aggregate of functions), the kernel and acquisition function (and their parameters) can be appropriately selected to guarantee convergence to the global optimum of the utility function.
Subsequently, \cref{subsection:globalconvergenceofDBAY} has shown that, if the agents are able to find the global optima of the aggregate of the utility functions in finite time, \DBay will converge to the global optimum of the objective function.
Combining these results proves the convergence of \DBay to the global optimum of the objective function for utility functions with known Lipschitz constants.

\section{Sample selection for continuous DCOPs}
\label{section:BayesianOptimization}

\structure{Which sample selection process to use?}
As mentioned in \cref{section:dcop_solvers}, in order to efficiently solve a continuous DCOP a sample-based solver is required.
To avoid discretization, a sample-based solver will be introduced such that the relation between the assignment and the utility will be taken into account within the sample selection process.

\structure{Overview of possibilities}
Within the literature several suitable methods exist for the sampling of continuous domains:
\begin{enumerate*}[label={\roman*)}]
    \item simulated annealing \cite{cerny1985}, which is a method that resembles the annealing of metals during cooling,
    \item genetic algorithms \cite{Goldberg1989}, which are a class of algorithms that resemble optimization through survival of the most \emph{fit} solutions, and
    \item Bayesian optimization \cite{Mockus1982}, which is a method to find the optimum of a function by using a probabilistic model of that function.
\end{enumerate*}
All methods mentioned above could potentially be used within a sample-based DCOP solver.
In this work \BO will be used as for certain probabilistic models convergence to the global optimum can be guaranteed as shown in \citeA[Theorem~6]{Vazquez2010}.
These results are not available for either simulated annealing or genetic algorithms.

Indeed, as noted by \citeA{Rios2013}, based on \citeA{Romeo1991}, although convergence results exist for simulated annealing no performance guarantees can be given for finite number of iterations.
Genetic algorithms, as summarized in \shortciteA{El-Mihoub2006}, are capable of converging to a region in which the global optimum exists; however they may require a large number of samples to converge to an optimum within this region \cite{de2005genetic}.
As highlighted by \citeA{Lee2018}, this optimum is not guaranteed to be the global optimum.

\structure{What is BO in detail?}
\BO is a method to find the global optimum of a function in a sample-efficient manner, i.e.\ it minimizes the number of required samples.
\BO consists of two elements: a \emph{probabilistic model} to approximate a (utility) function \mathcmd{\utilityFunction(\cdot)}, and an \emph{acquisition function} \mathcmd{\acquisitionFunctionBase(\cdot)} to optimally select a new sample \mathcmd{\observationInputIndexed}, where \mathcmd{s} denotes the sample index.
These two elements are discussed in more detail in the subsequent sections.
Based on samples of the function and the corresponding function values \functionOutputLongIndexed, the probabilistic model is used to estimate a mean function $\mathcmd{\mu}(\cdot)$ and the corresponding variance function $\mathcmd{\deviationFunctionBase^2}(\cdot)$.
Every input/output pair, \observationLongIndexed, is included in the ordered observations set \mathcmd{\observationSet_\mathcmd{S} = \{\mathcmd{O}_1, \dots, \mathcmd{O}_\mathcmd{S} \}}, where \mathcmd{S} is the number of observations.
The observations are (re)ordered after every new observation, such that for scalar arguments \mbox{$\mathcmd{x}_1 \leq \mathcmd{x}_2 \leq \dots \leq \mathcmd{x}_\mathcmd{S}$}.
The observations are used to update the probabilistic model, such that after every new observation the approximation is refined.
An overview of the Bayesian optimization algorithm is given in Algorithm \ref{algorithm:BO}.
\begin{algorithm}
    \SetKwInOut{Input}{Input}
    \SetKwInOut{Output}{Output}
    \Input{\mathcmd{\utilityFunction(\cdot)}, \mathcmd{\acquisitionFunctionBase(\cdot)}, \mathcmd{\bm{x}}}
    \Output{$\mathcmd{\mu}(\mathcmd{\bm{x}}, \mathcmd{\mathcal{O}})$, $\mathcmd{\deviationFunctionBase^2}(\mathcmd{\bm{x}}, \mathcmd{\mathcal{O}})$}
    \For{s = 1, 2, \dots}{
        \tcc{Select the next sample based on acquisition function}
        $\mathcmd{\observationInputIndexed} \coloneqq \argmax_{x} q(x | \mathcal{O}_{s-1})$\;
        \tcc{Sample the utility function}
        $y_s \coloneqq \mathcmd{f}(\mathcmd{\observationInputIndexed})$\;
        \tcc{Augment (and reorder) the observation set}
        $\mathcal{O}_{s} \coloneqq \mathcal{O}_{s-1} \cup \{ (\mathcmd{\observationInputIndexed}, y_s) \}$\;
        \tcc{Calculate the mean function and the variance function}
        $\mathcmd{\mu}(\mathcmd{\bm{x}}, \mathcmd{\mathcal{O}}) =
        \mathcmd{\mathds{E}}\left[ \mathcmd{f}(\mathcmd{\bm{x}}) | \mathcmd{\mathcal{O}} \right]
        $\;
        $\mathcmd{\deviationFunctionBase^2}(\mathcmd{\bm{x}}, \mathcmd{\mathcal{O}}) =
        \mathcmd{\mathds{E}}\left[
            \left(
                [\mathcmd{f}( \mathcmd{\bm{x}}) | \mathcmd{\mathcal{O}}] - \mathcmd{\mu}(\mathcmd{\bm{x}}, \mathcmd{\mathcal{O}})
            \right)^2
        \right]
        $\;
    }
    \caption{Bayesian optimization \cite{Mockus1982}}
    \label{algorithm:BO}
\end{algorithm}

\subsection{Probabilistic model}
\label{subsection:probabilistic_model}
The \GP{} is a widely used probabilistic model to represent acquired \mbox{knowledge} about a function.
More elaborate models exist, but these will often undo the computational benefit of the \GP model.
Using the \GP model, a function \mathcmd{\utilityFunction(\cdot)} is modeled based on a prior mean function \gpMeanFunctionLong and a kernel \kernelLong.
The kernel represents the cross-correlation between the two variables \mathcmd{x}, $\mathcmd{x}^\prime$.
The prior mean function and the kernel contain all (prior) knowledge of \mathcmd{\utilityFunction(\cdot)}.
Typically, the prior mean function is equal to the zero function ($\mathcmd{m}(\mathcmd{x}) = 0$ for all \mathcmd{x}) if no prior information about the function is available.
Therefore, the modelling of the function depends mostly on the choice of \mathcmd{\kappa}.

The \GP model is combined with the observations in order to construct the joint Gaussian distribution over the function.
From the joint Gaussian distribution, the posterior (distribution) can be found by using the Sherman-Morrison-Woodbury formula \cite{Sherman1950}:
\begin{align}
    \mathcmd{P}(\mathcmd{f}(\mathcmd{x}) | \mathcmd{\mathcal{O}}) &= \mathcmd{\mathcal{N}}\bracket{\meanFunction, \mathcmd{\varianceFunctionBase(\mathcmd{x} | \observationSet)} }, \label{eq:posterior}
\intertext{where}
    \meanFunction &= \meanFunctionExtended \label{eq:mean_function} \\
    \mathcmd{\varianceFunctionBase(\mathcmd{x} | \observationSet)} &= \varianceFunctionExtended \label{eq:variance_function}
\end{align}
and \mathcmd{\mathcal{N}} denotes the normal distribution,
$\mathcmd{\bm{K}}(\mathcmd{\mathcal{O}})$ is the Gramian matrix of the kernel, defined by \kernelGramianLong,
\kernelCrossCorrelationLong denotes the \emph{cross-correlation} vector between the observations and \mathcmd{x},
and \observationOutputsLong denotes the observation value vector.
The (posterior) mean and variance function of the probabilistic model are denoted as \meanFunctionDot and \mathcmd{\varianceFunctionBase(\cdot | \observationSet)}, respectively.
Note that the posterior distribution contains the estimate of the function based on both the prior knowledge and the observations.

A wide range of kernels for Gaussian processes exist in the literature and the interested reader is referred to the work of \citeA{Duvenaud2011} for a overview on constructing kernels.
An important kernel property is universal approximation, which indicates the ability to estimate every continuous function up to a required resolution given a sufficient number of observations.
A kernel that possesses this property is called a universal kernel.
In the work of \citeA{Micchelli2006a}, the conditions for a kernel to be universal in terms of properties of its features are given.
Additionally, several examples of kernels are given that possess this property.
The most commonly used universal kernel is the \SEKernel.
A general description of the kernel and its properties is given by \citeA{Vert2004}.
A drawback of the \SEKernel is that it can result in over-smooth approximations.
For this reason, the \Matern kernel \cite{Minasny2005} is often used, since it can trade off differentiability and smoothness.
In practice, the choice for a kernel depends on the properties of the problem that needs to be solved.
Additionally, depending on the kernel one or more parameters need to be set.
All kernels have parameters that can be used to adjust their properties, such as smoothness and scaling.
If information about \mathcmd{\utilityFunction(\cdot)} is available, this should be incorporated in the selection of the kernel and its parameters.
Typically, it is assumed that no information about \mathcmd{\utilityFunction(\cdot)} is available, and then, as noted by \citeA{Rasmussen2004}, the selection of the parameters is non-trivial.
For this reason, the selection of parameters is often treated as a separate optimization problem \cite{MacKay1992}.
It is commonly solved by using the maximum likelihood criterion for which \ARD \cite{Mackay1994} is a widely used algorithm.

\subsection{Acquisition function}
\label{subsection:acquisition_function}
The selection of the next sample is the result of the optimization of an acquisition function \mathcmd{\acquisitionFunctionBase(\cdot)}, defined by
\begin{align*}
    \mathcmd{\observationInputIndexed} &= \argmax_{\mathcmd{x}} \mathcmd{q}(\mathcmd{x} | \mathcmd{\mathcal{O}}).
\end{align*}
The acquisition function depends on the posterior distribution in \cref{eq:posterior} and thereby on all previous observations.
Two commonly used acquisition functions are the \PI function \cite{Kushner1964} and the \EI function \cite{Mockus1978}.
The probability of improvement function considers the \emph{probability} of finding an observation the value of which is larger than the maximum observed value.
The maximum observed value is defined as
\begin{align*}
    \observationMaxSampleLong.
\end{align*}
The corresponding input is defined as \observationMaxSampleInputLong.
As noted by \citeA{brochu2010tutorial}, the probability of improvement function focusses solely on exploitation of already observed samples.
In order to balance exploration of the search space and exploitation of the observations, the expected improvement function will be used in this work.
The expected improvement function chooses the sample based on the expected \emph{value} of the next observation.
The expected improvement function can be written in closed form in terms of the mean and the variance function of the probabilistic model as
\begin{align}
    \acquisitionFunction &=
    \begin{cases}
        \acquisitionFunctionNonZeroCumulativeLong
        +
        \acquisitionFunctionNonZeroDensityLong
        &\text{if }
        \mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} > 0\\
        0 &\text{if } \mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} = 0
    \end{cases}
    \label{eq:acquisition}
    \\
    \acquisitionHelper &= \acquisitionHelperExtended
\end{align}
where $\Phi(\cdot)$ is the Gaussian cumulative distribution function, $\phi(\cdot)$ is the Gaussian probability density function, and \mathcmd{\xi} is a design parameter.
The design parameter can be used to tradeoff exploration and exploitation.
As noted by \citeA{Lizotte2012}, even for values as low as $\mathcmd{\xi} = 0$ will not result in a solely exploiting sampling behavior.
The interested reader is referred to \citeA{brochu2010tutorial} for a comparison of the two acquisition functions and more details.

\section{Conclusions}
\label{section:conclusion}
In this work, the novel algorithm called \DBayLong has been introduced to solve a \DCOPLong with continuous domains.
Within \DBay, the continuous domains are sampled based on \BO.
This removes the need for discretization of the domains.
Compared to traditional DCOP solvers, which require discretization, it results in a reduction in the computational demands of the individual agents.
For utility functions with known Lipschitz constants, \DBay is proven to converge to the global optimum solution of the DCOP.

A sensor coordination problem has been used to evaluate the performance of \DBay.
The results show that \DBay outperforms a centralized approach based on the achieved utility as a function of the required number of samples.
This sample efficiency is a result of the application of \BO with \DBay.

In future work, \DBay will be extended towards dynamic DCOPs \cite{fioretto2019distributedreview} in which the agents need to optimize a dynamic problem at every time step.
An extension will increase the applicability of the developed algorithm to dynamic real-world problems in which tracking of targets is an important factor, such as multi-agent surveillance.
In a dynamic adaptation of the sensor coordination problem, the locations of the targets change over time based on the target properties, such as speed and turn radius.
Furthermore, an extension towards non-ideal communication will be investigated.
Within that extension, the (pseudo-tree) communication structure will become dynamically generated and the message passing between the agents will be asynchronous.
An additional extension would be to further decrease the computational complexity of the algorithm by removing the optimization of the acquisition function.
A possible direction would be to apply a a bound-based method such as the approach of \citeA{Kawaguchi2016}.
\section{Distributed Constraint Optimization Problems}
\label{section:DCOP}
A \DCOPLong is a problem in which an objective function needs to be optimized in a distributed manner through value assignments for all variables.
The objective function consists of the aggregate of utility functions, which define the utility value of the value assignments.
Agents are able to assign a value for variables that are allocated to them.
Furthermore, a variable can only be allocated to a single agent.
Typically, the number of variables is equal to the number of agents, i.e. every agent assigns a single variable.
The agents cooperate by sending messages to agents with whom they share a utility function.
A utility function is shared between agents if their variables are in the arguments of that function.
An important aspect of a DCOP is the definition of the domains of the variables.
A domain defines all possible values a variable can be assigned to.
In other words, the value assignments are restricted by the domains of the variables.

Following the notation of \citeA{fioretto2019distributedreview}, a continuous DCOP is defined by \mbox{\mathcmd{\DCOPTuple = \langle\mathcmd{\mathbf{A}}, \mathcmd{\mathbf{X}}, \mathcmd{\mathbf{D}}, \mathcmd{\mathbf{F}}, \mathcmd{\alpha}, \operator\rangle}} where
\begin{itemize}
    \item \mathcmd{\agentSet = \{ \agent_\mathcmd{1}, \dots, \agent_\mathcmd{M} \}} is the set of agents,
        where \mathcmd{M} is the number of agents.

    \item \variableSetLong is the set of variables, where $\mathcmd{N} \geq \mathcmd{M}$ is the number of variables.

    \item \mathcmd{\domainSet = \{ \domainSet_\mathcmd{1}, \dots, \domainSet_\mathcmd{N} \}} is the set of domains of all variables, where $\mathcmd{\domainSet_\mathcmd{i}} \subseteq \mathcmd{\mathds{R}}$ is the (continuous) domain associated with variable \mathcmd{\variable_\mathcmd{i}}.
        The search space of the DCOP is defined by all possible combinations of all values within the domains as \searchSpaceLong, where $\prod$ is the set Cartesian product.
        The search space of a set of variables ($\mathcmd{\mathbf{V}} \subseteq \mathcmd{\mathbf{X}}$) is defined as $\mathcmd{\mathbf{\Sigma}}_\mathcmd{\mathbf{V}} = \underset{\mathcmd{i} ~:~ \mathcmd{\variable_\mathcmd{i}} \in \mathcmd{\mathbf{V}}}{\prod} \mathcmd{\domainSet_\mathcmd{i}}$.

        An assignment denotes the projection of variables onto their domain as \assignmentCompleteLong.
        In other words, for all $\mathcmd{\variable_\mathcmd{i}} \in \mathcmd{\mathbf{X}}$ if \assignmentIndexed is defined, then $\assignmentIndexed \in \mathcmd{\domainSet_\mathcmd{i}}$.
        An assignment of a subset of variables is denoted by \assignmentPartialLong.

    \item \mathcmd{\functionSet = \{ \utilityFunction_\mathcmd{1}, \dots, \utilityFunction_\mathcmd{K} \}} is the set of utility functions,
        where \mathcmd{K} is the number of utility functions.
        The scope of \mathcmd{\utilityFunction_\mathcmd{n}} is denoted as \variableSetFunctionIndexedLong, where $\mathcmd{\variable_\mathcmd{i}} \in \variableSetFunctionIndexed$ when \mathcmd{\variable_\mathcmd{i}} is an argument of \mathcmd{\utilityFunction_\mathcmd{n}}.
        The optimum of a utility function is denoted by \functionOutputOptimalLong with input \functionInputOptimalLong, where $\mathcmd{\searchSpace_{\variableSetFunctionIndexed}}$ denotes the domain of the utility function.

    \item \mathcmd{\variableMapping :~ \mathcmd{\mathbf{X}} \rightarrow \mathcmd{\mathbf{A}}} is a mapping from variables to agents.
        The agent to which variable \mathcmd{\variable_\mathcmd{i}} is allocated is denoted as \variableMappingAgentIndexed.
        A common assumption is that the number of agents is equal to the number of variables, such that $\mathcmd{\agent_\agentIndex} = \variableMappingAgentIndexed$ for $\mathcmd{i} = 1, \dots, \mathcmd{N}$.
        Likewise, the set of agents associated with \mathcmd{\utilityFunction_\mathcmd{n}} is denoted by \functionScopeLongIndexed.

    \item $\operator$ is an operator that combines all utility functions into the objective function.
        Common options are the \emph{summation} operator $\bracketBig{\sum(\cdot)}$ and the \emph{maximum} operator $\big(\max(\cdot) \big)$.
        The objective function is defined by \globalUtilityLong.
        The optimal value assignment is denoted by \assignmentOptimalLong.

\end{itemize}

The relation between the variables and the utility functions is typically represented as a constraint graph.
In this representation the agents are defined as nodes and the edges implicitly represent the utility functions.
A constraint graph is often converted into a pseudo-tree \cite{Freuder1985} through a distributed depth-first-search procedure to partition the problem into independent subproblems based on the utility functions.
A pseudo-tree is a rooted spanning tree that introduces hierarchy to the agents where the subproblems are contained in separate \emph{branches}.
Additionally, all agents are assigned a single parent, which is an agent higher in the hierarchy.
The only exception is the agent on top of the hierarchy, which is called the root of the tree, this agent has no parent.
An agent can have multiple children and the agents without children correspond to the leaves of the tree.
In addition to parent/child relations, the pseudo-tree defines pseudo-parent/pseudo-child relations to indicate relations between agents over multiple hierarchy levels.
Typically, the pseudo-tree is used as a communication structure, where agents only communicate between parent and child.
In these cases, the pseudo relation allows for (indirect) interaction between pseudo-parent and pseudo-children.
A graphical example of the two DCOP representations is given in \cref{fig:DCOP_representations}.

\begin{figure}[h]
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
        \input{figures/DCOP_representation_constraint_graph.tikz}
        \caption{Constraint graph}
        \label{fig:dcop_constraint_graph}
    \end{subfigure}
    \qquad \qquad
    \begin{subfigure}[b]{0.3\textwidth}
        \input{figures/DCOP_representation_pseudo_tree.tikz}
        \caption{Pseudo-tree}
        \label{fig:dcop_pseudo_tree}
    \end{subfigure}
    \caption{
        Graphical example for two representations of the same DCOP.
        A node represents an agent $\mathcmd{\agent_\agentIndex}$ and the edges indicate the utility functions $\mathcmd{\utilityFunction_\mathcmd{n}}$ based on the variables $\mathcmd{\variable_\mathcmd{i}}$.
        In \cref{fig:dcop_constraint_graph} the nodes are unstructured, while in \cref{fig:dcop_pseudo_tree} the hierarchy is indicated through separate layers.
        The agent at the top ($\mathcmd{a}_\mathcmd{1}$) is the root of the tree, and the bottom agents ($\mathcmd{a}_2, \mathcmd{a}_4, \mathcmd{a}_5$) are the leaves.
        Note that the edge between $\mathcmd{a}_\mathcmd{1}$ and $\mathcmd{a}_4$ specifies a pseudo-parent/pseudo-child relation.
        }
        \label{fig:DCOP_representations}
\end{figure}

\section{Distributed Pseudo-tree Optimization Procedure (DPOP)}
\label{section:DPOP}

The DPOP algorithm \cite{Petcu2005a} solves a DCOP in three phases:
\begin{enumerate}
    \item[(2)] \emph{Bottom-up utility propagation}:
        From the leaves (nodes without children), the agents pass a utility message \utilityMessageLongDomain towards their parents.
        The message represents the projection of the (partial) global utility of the utility functions of the agents and its children with respect to all variables of the parent (and pseudo-parents).
        In words, it is a multidimensional array with one dimension for each variable within the message where every element represents an allocation of the variables, $\mathcmd{\mathfrak{U}}[\mathcmd{\mathbf{X}}] = \operator(\mathcmd{\mathbf{F}}[\mathcmd{\mathbf{X}}])$.
        \\
        \\
        A utility message \mathcmd{\utilityMessage_\agentIndex^\agentIndexOther} is send by agent \mathcmd{\agent_\agentIndex} to agent \mathcmd{\agent_\agentIndexOther} based on the utility function set and utility messages received by its children.
        \\
        \\
        Messages are combined through the join operation $\mathcmd{\mathfrak{U}}^\agentIndex_{1,2} = \mathcmd{\mathfrak{U}}^\agentIndex_{1} \join \mathcmd{\mathfrak{U}}^\agentIndex_{2}$.
        This results in the union of the scopes, i.e.\
        $\mathcmd{\variableMapping}(\mathcmd{\mathfrak{U}}^\agentIndex_{1,2}) = \mathcmd{\variableMapping}(\mathcmd{\mathfrak{U}}^\agentIndex_1) \cup \mathcmd{\variableMapping}(\mathcmd{\mathfrak{U}}^\agentIndex_2)$, and sums the value of the elements for all combinations of $\mathcmd{\mathfrak{U}}^\agentIndex_{1,2}[\mathcmd{\mathbf{X}}] = \operator(\mathcmd{\mathfrak{U}}^\agentIndex_{1}[\mathcmd{\mathbf{X}}], \mathcmd{\mathfrak{U}}^\agentIndex_{2}[\mathcmd{\mathbf{X}}])$.
        \\
        \\
        All child messages are combined as \utilityMessageChildrenIndexedLong.
        The local message holds the utility values for the utility functions as, \utilityMessageLocalIndexedLong.
        The (pseudo-)parents message, which hold the utility values of the functions shared by the agent and its (pseudo-)parents, are combined as \utilityMessageParentsIndexedLong.
        \\
        \\
        The utility messages of the children \mathcmd{\utilityMessage^\agentIndex_{\Child_{\agentIndex}}} and (pseudo-)parents \mathcmd{\utilityMessage^\agentIndexOther_{\Parent_{\agentIndex}}} are combined with the utility message \mathcmd{\utilityMessage_\agentIndex} to form the total utility message \mathcmd{\utilityMessage^\agentIndex} of agent \mathcmd{\agent_\agentIndex} as \utilityMessageAgentLongIndexed.
        \\
        \\
        The resulting message \mathcmd{\utilityMessage^\agentIndex} is optimized over the local variables of the agent \mathcmd{\agent_\agentIndex}.
        This operation is defined by a projection operator $\project$ and
        assigns the allocation of \variableMappingAgentIndexed that results in the maximal utility.
        The result is a message of lesser cardinality based on the number of variables in \variableMappingAgentIndexed.
        Formally, \mathcmd{\mathcmd{\utilityMessage^\agentIndex} \perp \variableMappingAgentIndexed = \max_{\variableMappingAgentIndexed} \mathcmd{\utilityMessage^\agentIndex}[\variableMappingAgentIndexed]}, and the scope \mathcmd{\mathcmd{\variableMapping}(\mathcmd{\utilityMessage^\agentIndex} \perp \variableMappingAgentIndexed) = \mathcmd{\variableMapping}(\mathcmd{\utilityMessage^\agentIndex}) \setminus \variableMappingAgentIndexed}.
        \\
        \\
        For example, when $\mathcmd{\utilityMessage^\agentIndex} = \bigl[
            \begin{smallmatrix}
            4 & 5 \\
            6 & 7
            \end{smallmatrix}
            \bigr]
            $ with
            $\mathcmd{\variableMapping}(\mathcmd{\utilityMessage^\agentIndex}) = \{ \mathcmd{x}_1, \mathcmd{x}_2 \}$,
            the projection is
            $\mathcmd{\utilityMessage^\agentIndex} \project \{\mathcmd{x}_1\} =
            \bigl[
            \begin{smallmatrix}
            6 & 7 \\
            \end{smallmatrix}
            \bigr]$ and $\mathcmd{\mathfrak{X}}(\mathcmd{x}_1 | \mathcmd{x}_2) = \left\{ \{ \mathcmd{x}_1 = 2 | \mathcmd{x}_2 = 1 \}, \{ \mathcmd{x}_1 = 2 | \mathcmd{x}_2 = 2 \} \right\}$.
        \\
        \\
        After the projection, the message is sent towards to its parent,
        $\mathcmd{\mathfrak{U}}_i^{\mathcmd{\mathrm{\bf{P}}}_i} = \mathcmd{\utilityMessage^\agentIndex} \project \variableMappingAgentIndexed$.
        After the root agent of the pseudo-tree (agent without parent) has finished this procedure, the value propagation phase is initiated.

    \item[(3)] \emph{Top-down value propagation}:
        Based on its child messages, the root agent chooses the optimal assignment of its local variables \variableMappingAgentIndexed as,
        $\mathcmd{\allocation^*}_\agentIndex = \argmax_{\mathcmd{\mathfrak{X}}_\agentIndex}(\mathcmd{\utilityMessage^\agentIndex}[\variableMappingAgentIndexed])$.
        The allocation of these values is sent to all the children of the agent.
        Based on these values the children allocate their own variables,
        $\mathcmd{\allocation^*}_\agentIndexOther = \argmax_{\mathcmd{\mathfrak{X}}_\agentIndexOther}(\mathcmd{\mathfrak{U}}^\agentIndexOther[\variableMappingAgentOtherIndexed] \perp \mathcmd{\allocation^*}_\agentIndex)$.
        Every agent repeats this process until the leaves of the tree are reached, completing the assignment $\mathcmd{\allocation^*}$ for all agents.
\end{enumerate}
\section{DCOP solvers}
\label{section:dcop_solvers}

\structure{Overview of DCOP, originate from CSP, therefore discrete}
The DCOP framework originates from an extension and generalization of Constraint Satisfaction Problems (CSPs) \cite{Tsang1993} towards distributed optimization.
A solution for a CSP is defined as the assignment of all variables from (finite) discrete domains such that all constraints are satisfied.
The CSP framework has been extended from centralized optimization to agent-based distributed optimization in the work of \citeA{Yokoo1998}.
Within the Distributed-CSP framework, the variables are allocated to agents and the agents coordinate the assignments of the variables among each other.

Additionally, CSP has been generalized into the Constraint Optimization Problem (COP) framework, where the constraints are replaced with utility functions.
These utility functions return a cost or reward based on the assignment of the variables.
Instead of constraint satisfaction, the goal of a COP is to find assignments that optimize an objective function.
Finally, the DCOP framework provides a unified framework that includes a large class of problems by combining the generalization and the extension.
A graphical overview of the relations between the problem frameworks can be seen in \cref{fig:COP_DCOP}.
\begin{figure}[h]
    \centering
    \input{figures/COP_DCOP.tikz}
    \caption{
        Graphical overview of the relations between the problem frameworks.
        Adapted from \citeA{fioretto2019distributedreview}.
        }
        \label{fig:COP_DCOP}
    \end{figure}

However in the DCOP framework, the definition of the domains remains discrete, which limits its application to problems with continuous domains.
In the current paper, the (discrete) DCOP framework will be further extended to overcome this restriction in order to include real-world problems such as cooperative search.

\structure{How to extend/adapt the solvers for continious domains?}
As noted by \citeA{Modi2005}, optimally solving a discrete DCOP is NP-hard with respect to the number of variables and the cardinality of their domains.
For this reason, complete (optimal) DCOP solvers are often not used in practice.
In the literature a diverse range of incomplete (near-optimal) DCOP solvers exist that trade off solution quality against computational requirements.
Such solvers perform well for problems with domains a with low cardinality, such as graph coloring problems \cite{Modi2005} and meeting scheduling problems \cite{Jennings1995}.

Current DCOP solvers require the discretization of the continuous domains in the problem definition phase.
One typically discretizes the entire domain using a grid-based approach that converts the continuous domains into discrete domains.
This process can arbitrarily increase the cardinality of the domains, thereby rendering the use of current DCOP solvers intractable.
This raises the question whether there are DCOP solvers that can be extended to solve continuous DCOPs without the need for discretization.
\clearpage
In the literature numerous solvers for DCOP with discrete domains have been proposed; for a detailed overview, the reader is referred to \citeA{Cerquides2014} and \citeA{Leite2014}.
In order to determine if DCOP solvers can be extended towards continuous domains, the taxonomy introduced by \citeA{Yeoh2008} is used.
The taxonomy defines three classes:
\begin{description}
    \item [Search-based solvers] perform a distributed search over the local search space of the agents.
    These solvers are based on centralized search techniques such as best-first and depth-first to \emph{reduce the search space} of the problem by exchanging messages between the agents.
    Examples are ADOPT \cite{Modi2005}, CoCoA \cite{Leeuwen2017}, AFB \cite{Gershman2009}, DSA \cite{Kirkpatrick1983}, and DBA \cite{Wittenburg2003}.
    \item [Inference-based solvers] communicate accumulated information among agents in order to \emph{reduce the problem size after every message} through dynamic programming methods.
    Well-known examples of this class of solvers are DPOP \cite{Petcu2005a}, the max-sum based algorithm \cite{Rogers2011}, and action GDL \cite{Vinyals2009}.
    \item [Sampling-based solvers] coordinate the sampling of the global search space guided by probabilistic measures.
    The probabilistic measures are calculated based on (all) preceding samples in order to \emph{balance exploration and exploitation} of the global search space.
    At the time of writing, two sample-based solvers are found in the literature: DUCT \cite{Ottens2018} and Distributed Gibbs \cite{Nguyen2019}.
\end{description}

In order to efficiently solve a continuous DCOP the relation between the value of the variables and the corresponding (global) utility values should to be taken into account explicitly.
Neither search-based nor inference-based solvers can be extended to take this property into account because both solver types only compare the utility of value assignments within a single iteration.
On the contrary, sampling-based solvers take previous iterations into account when selecting an additional sample.
This feature can be extended to take the relation between the samples into account during the sample selection process.
Which would allow for selecting a sample from a continuous domain directly.
Therefore, in this work, a \emph{novel sampling-based solver} will be introduced that removes the need for discretization entirely.

\section{The Distributed Bayesian algorithm}
\label{section:distributed_bayesian}

\structure{Small introduction}
In the previous sections, background information has been given about the DCOP framework and \BO.
In this section, the novel sample-based DCOP solver \DBayLong is presented.
This solver is capable of directly solving continuous DCOPs without discretization of the domains.
\DBay uses \BO as the probabilistic measure to optimize the sample selection.
The overall procedure is similar to state-of-the-art sample-based solvers, e.g. DUCT \cite{Ottens2018} and \Gibbs \cite{Nguyen2019}.

\structure{Sample-based solvers in general}
Sample-based solvers coordinate the sampling of the global search space guided by probabilistic measures in order to balance exploration and exploitation.
The general outline of sample-based solvers is as follows.
Based on a pseudo-tree representation of the DCOP, the variables and utility functions are allocated to the agents.
Afterwards, two consecutive phases are iteratively repeated until termination conditions are satisfied.
The first phase, the sampling phase, is top-down and starts from the root agent.
The root starts the sampling phase by selecting a sample for all its variables.
A sample can be viewed as a temporary assignment of a variable.
The sample is sent as a \mathcmd{\mathrm{\bf sample}} message to all the children of the agent.
Upon receiving this message, an agent samples its own variables and adds these samples to the \mathcmd{\mathrm{\bf sample}} message before sending it to its own children.
This process continues until the leaf agents are reached.

When the leaf agents are reached, the utility phase is initiated.
This second phase is bottom-up and starts from the leaf agents.
Based on the allocated functions, the agents calculate the utility based on the \mathcmd{\mathrm{\bf sample}} message and the assignments of their own variables.
This utility is encoded within a \mathcmd{\mathrm{\bf utility}} message and sent to the parent of the agent.
Upon receiving a \mathcmd{\mathrm{\bf utility}} message, an agent calculates the utility of its utility functions.
The resulting utility value is added to the \mathcmd{\mathrm{\bf utility}} message before sending it to its parent.
This phase finishes when the root agent receives a \mathcmd{\mathrm{\bf utility}} message from all its children.
This moment marks the end of a single iteration within the solver.
At this time, all agents have obtained the utility value associated with the sample of their variables.
This information is used to update the probabilistic measure and thereby the selection of the sample in the next iteration.

\structure{Sample-based solver alternatives}
While the overall procedure is similar for DUCT and \Gibbs, their main difference is in the method of selecting additional samples.
The probabilistic measure in DUCT calculates confidence bounds of the utility of the samples and selects samples to improve these bounds.
The \Gibbs algorithm selects samples based on joint probability distributions and keeps track of the differences between the utility values of the samples as termination criterion.
While \Gibbs is more memory efficient compared to DUCT, both algorithms have a runtime complexity that is linear in the cardinality of the largest domain \cite[Table 4]{fioretto2019distributedreview}.
Therefore, both \Gibbs and DUCT suffer from discretization of continuous domains and are not suitable for continuous DCOPs.

\structure{Consecutive sampling and utility phase result in a problem}
An additional disadvantage of both solvers is the non-determinism with respect to the utility value of a sample.
This is caused by the consecutive sampling and utility phases, since within an iteration all agents sample a single value from their local search space.
In other words, the same \mathcmd{\mathrm{\bf sample}} message can result in different \mathcmd{\mathrm{\bf utility}} messages when the children of an agent select different samples for their variables.

In order to remove the non-determinism, the sampling and utility phase in \DBay will be restricted to parent and children instead of the entire pseudo-tree.
To be more precise, when a child receives a \mathcmd{\mathrm{\bf sample}} message it will first iterate between its own children before sending a \mathcmd{\mathrm{\bf utility}} message to its parent.
This will guarantee that the utility value of the \mathcmd{\mathrm{\bf utility}} message will always be the same for the same \mathcmd{\mathrm{\bf sample}} message.

\structure{Formal algorithm description}
\DBay as described in \cref{algorithm:DBAY} (\cref{section:distributed_bayesian_algorithm}) involves four phases:
\begin{description}[align=left]
    \item[(1) Pseudo-tree construction]
        The agents create a pseudo-tree from the constraint graph of the DCOP, by performing a depth-first search traversal \cite{Awerbuch1985}.
        Afterwards, every agent \mathcmd{\agent_\agentIndex} knows its parent/children sets (\mathcmd{\Parent_{\mathcmd{i}}}/\mathcmd{\Child_{\mathcmd{i}}}) and pseudo-parents/pseudo-children sets (\mathcmd{\PseudoParent_{\mathcmd{i}}}/\mathcmd{\PseudoChild_{\mathcmd{i}}}), where $\mathcmd{\Parent_{\mathcmd{i}}}, \mathcmd{\PseudoParent_{\mathcmd{i}}}, \mathcmd{\Child_{\mathcmd{i}}}, \mathcmd{\PseudoChild_{\mathcmd{i}}} \subset \mathcmd{\mathbf{A}}$.
        The pseudo-tree is used as communication structure in which agents only communicate with agents with whom they share a parent/child relation.
        Note that (depending on the DCOP) none of the agents have a complete overview of the pseudo-tree.

    \item[(2) Allocation of utility functions]
        Similar to the allocation of variables, all utility functions in \mathcmd{\mathbf{F}} are exclusively allocated to the agents.
        Every agent \mathcmd{\agent_\agentIndex} constructs two separate function sets based on the variables of the agent and the variables of its (pseudo-)parents.
        Firstly, the utility function set \functionSetLocalLongIndexed, which only depends on the agent itself.
        Secondly, the shared utility function set, \functionSetSharedLongIndexed, which involves the agent and its (pseudo-)parents.
        These two function sets are combined as \functionSetLocalCombinedIndexedLong.

    \item[(3) Sample propagation]
        In this phase, every agent optimizes its local variables through the \BO method and the exchange of \mathcmd{\mathrm{\bf sample}} and \mathcmd{\mathrm{\bf utility}} messages.
        By doing so, the assignments of the local variables will converge to the global optimum of the objective function as will be shown in \cref{subsection:globalconvergenceofDBAY}.
        The local variables of \mathcmd{\agent_\agentIndex} are defined as \localVariablesIndexedLong.
        The variables are optimized based on the aggregate utility of all utility functions in set \functionSetLocalCombinedIndexed and the utility functions of its children ($\mathcmd{\utilityFunction_\mathcmd{n}} \in \functionSetLocalCombinedChild$ for all $\mathcmd{a}_\agentIndexChild \in \mathcmd{\Child_{\mathcmd{i}}}$).
        Consequently, the aggregate utility values obtained by the root agent hold the utility values of the objective function.

        Since a sample from (pseudo-)parents is required to calculate the utility of the functions in \functionSetSharedIndexed, every agent \mathcmd{\agent_\agentIndex} waits for a \mathcmd{\mathrm{\bf sample}} message.
        The phase is therefore initiated by a \mathcmd{\mathrm{\bf sample}} message from the root agent.
        The phase finishes when a convergence threshold is reached by the root agent.
        Upon receiving \mathcmd{\mathrm{\bf sample}} message \mathcmd{\sampleMessage_\agentIndexParent} from its parent \mathcmd{\agent_\agentIndexParent}, agent \mathcmd{\agent_\agentIndex} optimizes its own variables accordingly by sampling all functions in set \functionSetLocalCombinedIndexed.
        The samples are selected through the optimization of an acquisition function.
        Note that the acquisition function is based a kernel and on all preceding samples.

        If the agent is a leaf agent, the agent can optimize its variables without considering the impact of its assignments on other agents.
        However, when the agent has children, the agent needs to send its sample and wait for a \mathcmd{\mathrm{\bf utility}} message to be able to know its effect on the functions in the sets \functionSetLocalCombinedChild  for all $\mathcmd{a}_\agentIndexChild \in \mathcmd{\Child_{\mathcmd{i}}}$.
        Therefore, the agent augments the \mathcmd{\mathrm{\bf sample}} message of its parent with its own sample as \sampleMessageIndexedOtherLong and sends the message to all its children.
        The agent then waits until it has received all \mathcmd{\mathrm{\bf utility}} messages from its children.
        Only then the agent is able to compute the aggregate utility and return a \mathcmd{\mathrm{\bf utility}} message to its own parent.
        Note that the aggregate utility represents the optimal utility for the sample of the agent and all its \mbox{(pseudo-)children}.
        A \mathcmd{\mathrm{\bf utility}} message is defined as \optimumMessageToParentLong, where \mbox{\optimumMessageLocalLong} and \optimumMessageIndexedChildren define the utility and the aggregated child utility, respectively.
        A graphical overview of the \emph{sample phase} is shown in \cref{fig:B-DPOP_sample_phase}.

    \item[(4) Assignment propagation]
        The final phase is the assignment propagation phase, in which the root agent $\mathcmd{a}_\mathcmd{1}$ sends the final assignment of all its variables to its children as a \mathcmd{\mathrm{\bf final}} message \sampleMessageBestInitialRootLong.
        Based on these assignments the children can assign their own variables to the value corresponding to the optimal utility value.
        Afterwards, every agent adds its assignments to the message as \sampleMessageBestInitialIndexedLong.
        After the leaf agents have received a \mathcmd{\mathrm{\bf final}} message, all agents have completed their local assignments \assignmentPartialAgentOptimal.
        Note that typically no agent has information of the complete assignment, \mbox{$\assignmentPartialOptimal = \{ \assignmentPartialAgentOptimal ~:~ \mathcmd{i} = 1, \dots, \mathcmd{M} \}$}.
\end{description}
\begin{figure}[ht]
    \centering
        \input{figures/Algorithm_sample_phase_example.tikz}
    \caption{
        Graphical overview of the sample phase of \DBay.
        Agents are indicated by circles labeled with an agent index, and utility functions are shown as black lines.
        Starting from the root $\mathcmd{a}_\mathcmd{1}$ (top-left), a \mathcmd{\mathrm{\bf sample}} message $\mathcmd{\mathfrak{S}}_\mathcmd{1}$ is sent to its children ($\mathcmd{a}_2, \mathcmd{a}_3$).
        Subsequently, agent $\mathcmd{a}_3$ will send a \mathcmd{\mathrm{\bf sample}} message $\mathcmd{\mathfrak{S}}_3$ to its children.
        After iterating between its children and calculating its local utility, agent $\mathcmd{a}_3$ combines all local utilities and sends a \mathcmd{\mathrm{\bf utility}} message $\mathcmd{\mathfrak{U}}^\mathcmd{1}_3$ to its parent.
        This process is repeated when the root $\mathcmd{a}_\mathcmd{1}$ sends another  \mathcmd{\mathrm{\bf sample}} message and finishes when $\mathcmd{a}_\mathcmd{1}$ reaches a convergence threshold.
        Note that the interactions between $\mathcmd{a}_\mathcmd{1}$ and $\mathcmd{a}_2$ are not illustrated.
    }
    \label{fig:B-DPOP_sample_phase}
\end{figure}

\clearpage
\structure{Bridge to analysis section}
In the next section the convergence of \DBay to the global optimum of a continuous DCOP is analysed.
\DBay utilizes \BO for the sample selection within the sample propagation phase.
For that reason, the performance of \DBay depends highly on the properties of the \BO method.
As mentioned in \cref{section:BayesianOptimization}, \BO consists of the combination of a kernel and an acquisition function.
Therefore, the analysis is focussed on the selection of the kernel, the acquisition function, and their parameters.

\section{Introduction}
\label{section:introduction}

\structure{CONTEXT: Real world problems such as UAV: distributed, communication, limited processing power and memory}
Many real-world problems can be modeled as multi-agent problems in which agents need to assign values to their local variables to optimize a global objective based on utility.
Examples include scheduling \cite{Sato2015}, mobile sensor coordination \cite{Zivan2015}, hierarchical task network mapping \cite{Sultanik2007}, and cooperative search \cite{Acevedo2013}.
Even though numerous algorithms exist that solve these problems, applying them in practice is often problematic, as complications arise from limitations in communication, computation, and/or memory.

\structure{These problems are well suited to be modeled within the DCOP framework}
The \DCOPLong framework is well suited to model the above-mentioned problems (as detailed in \citeA{Meisels2008,Modi2005,Petcu2005b,Gershman2009,Yeoh2012}).
Within the DCOP framework a problem is defined based on variables and on utility functions that are aggregated into an objective function.
Additionally, agents assign values to all the variables that are allocated to them.
Agents are considered neighbors if their variables are arguments of the same utility function.
Neighbors then cooperatively optimize their utility functions through the exchange of messages.
Within the DCOP framework the variables are constrained by their domains.
In other words, a domain defines all possible value assignments of a variable.
This explicit definition of the domains of the variables is a major benefit for real-world problems that are (input) constrained.
Many solvers developed for DCOP assume that these domains are finite and discrete, while real-world problems are typically characterized by finite continuous domains.
A common approach for DCOP solvers is to use equidistant discretization, such as using a grid overlay to define all possible positions of an agent in an area.
When using equidistant sampling to discretize a continuous domain, the quality of the solution will depend on the distance between the values, where a smaller distance will allow for a better solution.
However, when continuous domains are discretized their cardinality will grow.
The increase in cardinality will exponentially increase the search space.
From the overview articles of \citeA{Leite2014} and \citeA{fioretto2019distributedreview}, it is clear that the cardinality of the domains is a major restriction to current DCOP solvers.
Therefore, discretization can cause solving continuous DCOPs to be intractable for current (discrete) DCOP solvers despite a small number of variables.

\structure{Solution is to take relation into account}
The underlying reason for the increase in problem size is that current DCOP solvers (implicitly) consider all values within a domain as unrelated to each other.
Because of this assumption it is not possible to efficiently sample the search space.
In problems with continuous domains this assumption can be removed since the utility of values that are close is often similar.
By explicitly taking this relation into account, a DCOP can be solved using efficient optimization methods.

\structure{Of possible approaches to take this relation into account, BO = best}
Various types of optimization methods exist in the literature that take advantage of the relation between the value and the utility.
Examples include simulated annealing \cite{cerny1985}, genetic algorithms \cite{Goldberg1989}, and \BO \cite{Mockus1989}.
In this work, \BO will be used as it focusses on efficient sampling during optimization, thereby requiring relatively few samples to closely approach the optimum.

\structure{Contributions of work}
Overall, the contributions of this article are threefold.
Firstly, we introduce an efficient algorithm that uses methods found in \BO to solve DCOPs with continuous domains called \DBayLong.
Secondly, we provide a theoretical proof of the convergence of the proposed algorithm called \DBay to the global optimum of the DCOP for utility functions with known Lipschitz constants.
Lastly, simulation results are given for a sensor coordination problem to compare the sample efficiency of the solver to a centralized approach based on equidistant discretization of the continuous domains.

\structure{Summary of article}
The remainder of this paper is organized as follows.
Firstly, in \cref{section:dcop_solvers} relevant literature regarding DCOP solvers are discussed.
Background information about the DCOP framework and the Bayesian optimization algorithm is provided in Sections \ref{section:DCOP} and \ref{section:BayesianOptimization}, respectively.
Afterwards, we present the novel sampling-based DCOP solver called \DBay in \cref{section:distributed_bayesian}.
The theoretical properties of \DBay are analyzed in \cref{section:analysis}.
Evaluation of \DBay for a sensor coordination problem is included in \cref{section:simulation_results}.
Finally, \cref{section:conclusion} summarizes the results and defines future work.
\section{Simulation Results}
\label{section:simulation_results}
\newcommand{\mathcmd{\beta}}{\mathcmd{\beta}}
\newcommand{\mathcmd{l}}{\mathcmd{l}}
\newcommand{\mathcmd{\omega}}{\mathcmd{\omega}}

In this section, the performance of \DBay is empirically evaluated using two metrics.
The first is the achieved relative utility as a function of the number of samples.
The number of samples is used as the threshold for the agents in the sample phase.
The relative utility is defined as the achieved utility divided by the optimal utility.
The relative utility allows for the comparison of the results over various randomly generated problems.
The second metric focuses on the sample efficiency of the algorithm.
This metric is important to consider if the evaluation of the utility functions is computationally expensive.
This metric compares the number of samples required by \DBay and the number of samples required by the centralized approach to achieve equal relative utility.
In the centralized approach, the continuous domains are equidistantly discretized.
The optimum of all possible combinations of these discrete domains is then taken as the achieved centralized utility.

\subsection{Sensor Coordination Problem}
The performance of \DBay is evaluated based on the sensor coordination problem.
The sensor coordination problem is an optimization problem in which every agent needs to orient its sensor in order to observe targets as accurately as possible.
This problem is modeled within the DCOP framework as a distributed problem.
A real-world analogue would be the optimization of the orientation of multiple cameras based on image recognition.
The image recognition process can require significant computational effort depending on the image quality and the type of target.
This makes it computationally intensive to check every orientation of the camera image, especially for a centralized approach.

Within the sensor coordination problem, all sensors are identical in terms of their sensor range $\mathcmd{l}$ and angle of view $\mathcmd{\beta}$.
These properties, combined with the position of the sensor, determine the observation domain of the sensor.
This domain defines all locations that could be observed by the sensor.
The orientation $\mathcmd{\omega}_\mathcmd{i}$ of the sensor \mathcmd{i} changes the observed area within the observation domain.
A target is detected when it is located within this area.
For every detected target, positive utility is allocated to the agent responsible for the sensor.
The maximum utility is allocated when the sensor is oriented directly at the target.
The utility value decreases linearly towards the edges of the observation area.
The optimal utility is determined by the optimal solution of a centralized approach with \mbox{\num{720} samples} for every domain.
The parameters of the problem are the number of targets, the number of sensors, the sensor range \mathcmd{l}, and the angle of view of the sensors \mathcmd{\beta}.
The sensors are arranged in an equally distanced rectangular grid.
The distance between the sensors is such that the combined observation domains of all sensors is maximized without allowing unobservable areas between the sensors.
While the locations of the targets $t$ are randomly distributed within the combined observation domains of the sensors.
In the simulations, all problems are generated with \num{6} sensors, \num{12} targets, and with identical sensor properties where the sensor range is set to $\mathcmd{l} = 1$ and the angle of view is set to $\mathcmd{\beta} = \SI{36}{\degree}$.
A graphical example of the simulated sensor coordination problems can be seen in \cref{fig:sensor_coordination_problem_example}.

\begin{figure}[h]
	\centering
		\input{figures/Problem_overview.tikz}
	\caption{
        Graphical example of a sensor coordination problem with \num{6} sensors and \num{12} targets.
        The sensors $\mathcmd{\agent_\agentIndex}$ are arranged in a grid.
        The distance between the sensors is based on the sensor range \mathcmd{l}.
        The sensor range is indicated by a dotted circle centered around the position of the sensor.
        The observed area of the sensors is shown as shaded areas, and is based on the angle of view \mathcmd{\beta} and the orientation $\mathcmd{\omega}_\agentIndex$.
        The targets are denoted as annotated black circles.
    }
    \label{fig:sensor_coordination_problem_example}
\end{figure}

The sensor coordination problem is described within the continuous DCOP framework as follows:
\begin{itemize}
    \item \mathcmd{\agentSet = \{ \agent_\mathcmd{1}, \dots, \agent_\mathcmd{M} \}} is the set of agents,
        where \mathcmd{M} is the number of agents.
        The position of an agent is denoted as $p_\agentIndex \in \mathcmd{\mathds{R}}^2$.

    \item $\mathcmd{\mathbf{X}} = \{ \mathcmd{\omega}_1, \dots, \mathcmd{\omega}_\mathcmd{N} \}$ is the set of sensor orientations, where $\mathcmd{N} = \mathcmd{M}$.

    \item \mathcmd{\domainSet = \{ \domainSet_\mathcmd{1}, \dots, \domainSet_\mathcmd{N} \}}, where $\mathcmd{\mathbf{D}}_\agentIndex = (\SI{-180}{\degree}, \SI{180}{\degree})$ for all $\mathcmd{i} = 1, \dots, \mathcmd{N}$ indicating all possible orientations of the sensor.

    \item $\mathcmd{\mathbf{F}} = \{f_{\mathcmd{n}}\}_{\mathcmd{n}=1}^{\mathcmd{K}}$ is the set of utility functions associated with the observation of the targets.
    The number of targets is denoted by $T \in \mathds{N}$.
    A target is located at position $t_\mathcmd{n} \in \mathcmd{\mathds{R}}^2$.
    The utility functions of the targets are described as,
    $f_{\mathcmd{n}} = \max_{\mathcmd{i} = 1, \dots, \mathcmd{N}}\left( f_{\mathcmd{n},\mathcmd{i}} \right)$ for $\mathcmd{n} = 1, \dots, \mathcmd{K}$, where
    \begin{align*}
        f_{\mathcmd{n},\agentIndex} =
        \begin{cases}
            1 - | \mathcmd{\omega}_\mathcmd{i} - \angle \overrightarrow{p_{\agentIndex}t_{\mathcmd{n}}} | / \mathcmd{\beta}
            &\mbox{if } \| \overrightarrow{p_{\agentIndex}t_{\mathcmd{n}}} \| \leq \mathcmd{l} \mbox{ and } | \mathcmd{\omega}_\mathcmd{i} - \angle \overrightarrow{p_{\agentIndex}t_{\mathcmd{n}}} | \leq \mathcmd{\beta}
            \\
            0 &\mbox{otherwise}
        \end{cases}
    \end{align*}
    and $\overrightarrow{p_{\agentIndex}t_{\mathcmd{n}}}$ denotes the vector between the location of the target $t_{\mathcmd{n}}$ and the position of the agent \mathcmd{\agent_\agentIndex}.

    \item $\mathcmd{\alpha}(\mathcmd{\omega}_\mathcmd{i}) = \mathcmd{\agent_\agentIndex}$ for $\mathcmd{i} = 1, \dots, \mathcmd{N}$ allocating a single sensor to every agent.

    \item $\operator = \sum(\cdot)$, resulting in the goal function $G(\cdot) = \sum_{f_\mathcmd{n} \in F} f_\mathcmd{n}(\cdot)$.

\end{itemize}

\subsection{Results}
The performance results of \DBay compared to the centralized approach are presented in \cref{fig:simulation_results}.
This figure shows the average results for 30 randomly generated problems for 6 sensors and 12 targets.
The results show an increase in the achieved relative utility of \DBay compared to the centralized approach based on the number of samples.
The difference in achieved utility can be explained by investigating the sampling strategies.
The centralized approach samples the sensor orientations equidistantly.
Therefore, as the number of samples is increased, the resolution of the samples decreases uniformly for the centralized approach.
While \DBay samples dynamically to balance exploration and exploitation based on all previously acquired observations.
Consequently, \DBay will initially focus on exploration and eventually focus on exploitation.
This behavior is clearly visible in \cref{fig:simulation_results_utility} in the range between \num{3} and \num{10} samples.
Within this range \DBay samples the sensor orientations equidistantly focussing on exploration.
This sampling behavior is identical to the centralized approach, which can be seen in the similarity in achieved utility.
At \num{11} number of samples the achieved utility of \DBay increases substantially.
This can be explained based on the angle of view of \SI{36}{\degree} of the sensors during the experiments.
At \num{10} samples the entire observation domain of a sensor was observed.
This enabled the switch to exploitation of the observations during the sample generation.
It clearly shows the advantage of the dynamic sampling of \DBay over equidistantly sampling.
The advantage is even more prominent when comparing the required number of samples by the centralized approach to achieve equal utility to \DBay, as shown in \cref{fig:simulation_results_samples}.

\newcommand{7cm}{7cm}
\begin{figure}[h]
    \centering
    \begin{subfigure}[t]{0.5\textwidth}
		\centering
		\includegraphics[height=7cm]{figures/simulation_results_utility.pdf}
        \caption{Achieved relative utility}
        \label{fig:simulation_results_utility}
	\end{subfigure}
	\begin{subfigure}[t]{0.45\textwidth}
		\centering
		\includegraphics[height=7cm]{figures/simulation_results_samples.pdf}
        \caption{Required number of samples}
        \label{fig:simulation_results_samples}
	\end{subfigure}
    \caption{
        Simulation results for randomly generated sensor coordination problems with \num{6} sensors and \num{12} targets.
        The figures show the average result of \num{30} randomly generated problems.
        \cref{fig:simulation_results_utility} shows the achieved utility of \DBay and the centralized approach relative to the optimum.
        Note that the achieved utility for both algorithms is equal for \num{3} samples, since the first \num{3} samples are equidistantly spaced for \DBay.
        \cref{fig:simulation_results_samples} shows the number of samples required by the centralized approach to achieve the same utility as \DBay.
    }
    \label{fig:simulation_results}
\end{figure}

\subsection{Convergence of \DBay based on global optima of utility functions}
\label{subsection:globalconvergenceofDBAY}
\structure{We go by the assumption that the agents are able to find the optimal assignments}
As shown in \cref{subsection:convergenceofBO}, all agents are able to find the global optimum of the aggregate utility of their utility functions and the utility functions of their children through Bayesian optimization.
In this section it will be shown that given the global optima of the utility functions, \DBay will find the global optimum of the objective function.

\structure{Leaf agents will produce optimal assignments and utility values given a sample}
Within the sample phase of \DBay, none of the agents can optimize their variables without interaction with other agents.
The interaction involves the sending of (top-down) \mathcmd{\mathrm{\bf sample}} messages and (bottom-up) \mathcmd{\mathrm{\bf utility}} messages.
Therefore, for the leaf agents, the optimization depends on their utility functions and the \mathcmd{\mathrm{\bf sample}} message of their parent, \mathcmd{\sampleMessage_\agentIndexParent}, as
\begin{align}
    \label{eq:optimal_leaf_agent_assignment}
    \hat{\mathcmd{\sigma}}_{\localVariablesIndexed}
    =
    \displaystyle\argmin_{\mathcmd{\sigma} \in \mathcmd{\mathbf{\Sigma}}_{\localVariablesIndexed}}
    \displaystyle\operator
    \Big(
        \optimumMessageLocal
    \Big)
    =
    \displaystyle\argmin_{\mathcmd{\sigma} \in \mathcmd{\mathbf{\Sigma}}_{\localVariablesIndexed}}
    \displaystyle\operator
    \left(
        \displaystyle\operator_{\mathcmd{\utilityFunction_\mathcmd{n}} \in \functionSetLocalCombinedIndexed}
        \Big(
            \mathcmd{\utilityFunction_\mathcmd{n}}(\assignmentPartial
            ~|~
            \mathcmd{\sampleMessage_\agentIndexParent}
        \Big)
    \right)
    \qquad \forall \mathcmd{\agent_\agentIndex} ~:~ \mathcmd{\Child_{\mathcmd{i}}} = \emptyset.
\end{align}

When the kernel and acquisition function are selected as detailed in \cref{subsection:convergenceofBO}, the assignment $\hat{\mathcmd{\sigma}}_{\localVariablesIndexed}$ is optimal with respect to the assignments of the (pseudo-)parents of agent \mathcmd{\agent_\agentIndex}, since $\mathcmd{\sampleMessage_\agentIndexParent} = \hat{\mathcmd{\sigma}}_{\mathcmd{\Parent_{\mathcmd{i}}}} \cup \hat{\mathcmd{\sigma}}_{\mathcmd{\PseudoParent_{\mathcmd{i}}}}$.
Consequently, the optimal assignment results in the optimal value for the \mathcmd{\mathrm{\bf utility}} message \optimumMessageToParent given the \mathcmd{\mathrm{\bf sample}} message \mathcmd{\sampleMessage_\agentIndexParent}.

\structure{All other agents therefore receive optimal utility values}
This optimal value is sent as a \mathcmd{\mathrm{\bf utility}} message to the parents of the leaf agents and results in the following assignment for the other agents:
\newcommand{\optimumMessageBestOptimal}{\mathcmd{
  \hat{\mathcmd{\mathfrak{U}}}^*
}}
\begin{align}
    \label{eq:optimal_agent_assignment}
    \hat{\mathcmd{\sigma}}_{\localVariablesIndexed}
    =
    \displaystyle\argmin_{\mathcmd{\sigma} \in \mathcmd{\mathbf{\Sigma}}_{\localVariablesIndexed}}
    \displaystyle\operator
    \Big(
        \optimumMessageLocal,
        \optimumMessageBestChild
    \Big)
    =
    \displaystyle\argmin_{\mathcmd{\sigma} \in \mathcmd{\mathbf{\Sigma}}_{\localVariablesIndexed}}
    \displaystyle\operator
    \left(
        \displaystyle\operator_{\mathcmd{\utilityFunction_\mathcmd{n}} \in \functionSetLocalCombinedIndexed}
        \Big(
            \mathcmd{\utilityFunction_\mathcmd{n}}(\assignmentPartial
            ~|~
            \mathcmd{\sampleMessage_\agentIndexParent}
        \Big)
        ,
        \optimumMessageBest_\agentIndex
    \right)
    \qquad \forall \mathcmd{\agent_\agentIndex} ~:~ \mathcmd{\Child_{\mathcmd{i}}} \neq \emptyset.
\end{align}
The aggregated \mathcmd{\mathrm{\bf utility}} message $\optimumMessageBest_\agentIndex = \displaystyle\operator_{\mathcmd{a}_\agentIndexChild \in \mathcmd{\Child_{\mathcmd{i}}}} \left( \optimumMessageFromChild \right)$ is the aggregate of the optimal \mathcmd{\mathrm{\bf utility}} messages of all children given an assignment of \mathcmd{\agent_\agentIndex}.
Therefore, agent \mathcmd{\agent_\agentIndex} is able to calculate the optimal assignment with respect to its parent \mathcmd{\mathrm{\bf sample}} message.

\structure{Finally, the root agent optimizes based on all aggregate \mathcmd{\mathrm{\bf utility}} messages}
This process is repeated until the root agent $\mathcmd{a}_\mathcmd{1}$ has received all \mathcmd{\mathrm{\bf utility}} messages from its children.
Since the root agent does not have any (pseudo-)parents, \cref{eq:optimal_agent_assignment} can be rewritten as
\begin{align}
    \hat{\mathcmd{\sigma}}_{\mathcmd{\mathcmd{\mathbf{X}}}_\mathcmd{1}}
    =
    \displaystyle\argmin_{\mathcmd{\sigma} \in \mathcmd{\mathbf{\Sigma}}_{\mathcmd{\mathcmd{\mathbf{X}}}_\mathcmd{1}}}
    \displaystyle\operator
    \Big(
        \mathcmd{\mathfrak{U}}_\mathcmd{1},
        \optimumMessageBest_\mathcmd{1}
    \Big)
    =
    \displaystyle\argmin_{\mathcmd{\sigma} \in \mathcmd{\mathbf{\Sigma}}_{\mathcmd{\mathcmd{\mathbf{X}}}_\mathcmd{1}}}
    \displaystyle\operator
    \left(
        \displaystyle\operator_{\mathcmd{\utilityFunction_\mathcmd{n}} \in \functionSetLocalCombinedRoot}
        \Big(
            \mathcmd{\utilityFunction_\mathcmd{n}}(\assignmentPartial)
        \Big)
        ,
        \optimumMessageBest_\mathcmd{1}
    \right)
    =
    \mathcmd{\sigma}^*_{\mathcmd{\mathcmd{\mathbf{X}}}_\mathcmd{1}}.
\end{align}
Note that $\optimumMessageBest_\mathcmd{1}$ holds the aggregate utility value of all other agents based on the sample of the root agent.
For that reason, if the root agent finds the optimal assignment $\hat{\mathcmd{\sigma}}_{\mathcmd{\mathcmd{\mathbf{X}}}_\mathcmd{1}}$ it is equal to the optimum of the objective function $\mathcmd{\sigma}^*_{\mathcmd{\mathcmd{\mathbf{X}}}_\mathcmd{1}}$.

After the root agent has found the optimal assignment of its variables it starts the final phase of \DBay.
In this phase the root agent sends the optimal assignment in a \mathcmd{\mathrm{\bf final}} message to its children, $\mathcmd{\hat{\sampleMessage}_\mathcmd{1}} = \{ \mathcmd{\sigma}^*_{\mathcmd{\mathcmd{\mathbf{X}}}_\mathcmd{1}} \}$.
Based on that optimal sample all agents are able to determine their optimal assignments, as shown in \cref{eq:optimal_agent_assignment}, and append their optimal assignment to the final message before sending the final message to their children, i.e. $\mathcmd{\hat{\sampleMessage}_\mathcmd{i}} = \{ \mathcmd{\sigma}^*_{\localVariablesIndexed} \} \cup \mathcmd{\hat{\sampleMessage}_\agentIndexParent}$.
This process is repeated until the leaf agents are reached and all agents have assigned the global optimal values to their variables.

\subsection{Convergence of \BO based on Lipschitz continuous functions}
\label{subsection:convergenceofBO}

\structure{Local optimization through dense sampling of the domains}
As shown by \citeA{torn1989global}, the convergence to the global optimum of a function by \BO can only be guaranteed through \emph{dense} sampling of the domain of the function.
For this reason, within the \BO method, the acquisition function will need to produce dense samples.
In the work of \citeA[Theorem~6]{Vazquez2010}, the \EI acquisition function, given by \cref{eq:acquisition}, is proven to produce dense observations within its search region.
The search region is defined in Definition \ref{def:searchRegion}.
\begin{defn}
    \label{def:searchRegion}
    The search region of the expected improvement acquisition function $\mathcmd{q}(\cdot, \mathcmd{\xi} | \mathcmd{\mathcal{O}})$ (based on \mathcmd{\utilityFunction(\cdot)} and \mathcmd{\mathcal{O}}) is defined by
    \begin{align*}
        \mathcmd{\searchRegion = \{ \mathcmd{x} \in \domainOperator(\mathcmd{f}) ~:~ \acquisitionFunctionBase(\mathcmd{x}, \acquisitionFunctionParameter | \observationSet) > 0 \}}.
    \end{align*}
\end{defn}

\structure{how to achieve dense sampling? -> how to ensure $S$ includes $x^*$?}
As a consequence of the dense sample generation property, \mathcmd{\mathcal{S}} will converge to an empty set when the number of samples goes to infinity.
Therefore, since the next sample is chosen from the search region ($\mathcmd{\observationInputIndexed} \in \mathcmd{\mathcal{S}}$), the global optimum will be sampled for $\mathcmd{s} \rightarrow \infty$ if $\mathcmd{\mathcmd{x}^*} \in \mathcmd{\mathcal{S}}$.
In other words, if the optimum inclusion ($\mathcmd{\mathcmd{x}^*} \in \mathcmd{\mathcal{S}}$) property holds, then global convergence is guaranteed.

\structure{Based on upper bound function the set $U$ includes $x^*$}
In order to find the conditions for which the optimum inclusion holds, the upper bound region is introduced.
The upper bound region set \mathcmd{\mathcal{U}} (Definition \ref{def:upperBoundRegion}) is based on the upper bound function \upperBoundFunctionDot (Definition \ref{def:upper_bound_function}).
Note that by definition, $\upperBoundFunctionBase(\mathcmd{x}) \geq \mathcmd{f}(\mathcmd{x})$ for all $\mathcmd{x}$.
As shown in Lemma \ref{lemma:optimumInU}, this region is guaranteed to include the global optimum if the optimum has not already been observed.
A graphical example of \mathcmd{\mathcal{U}}, \mathcmd{\mathcal{S}}, and \upperBoundFunctionDot can be seen in \cref{fig:LBEI}.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.85\textwidth, trim={0px 0px 0px 0px}, clip]{figures/LBEI.pdf}
    \caption{
        Graphical overview of the sets \mathcmd{\mathcal{U}}, \mathcmd{\mathcal{S}}, and the upper bound function \upperBoundFunctionDot based on the observations \mathcmd{\mathcal{O}}.
        Based on the observations (red circles), the corresponding upper bound function \upperBoundFunctionBase, the posterior mean function \mathcmd{\mu} and the the upper bound region \mathcmd{\mathcal{U}} are shown.
        }
        \label{fig:LBEI}
\end{figure}

\begin{defn}
    \label{def:upper_bound_function}
    The upper bound function \upperBoundFunction of (a Lipschitz continuous) function \mathcmd{\utilityFunction(\cdot)} over all observations in \mathcmd{\mathcal{O}} is defined by
    \begin{align*}
        \upperBoundFunctionLongIndex \qquad \forall \mathcmd{x} \in \domainOperator(\mathcmd{f}).
    \end{align*}
\end{defn}
\begin{defn}
    \label{def:upperBoundRegion}
    The upper bound region of \mathcmd{\utilityFunction(\cdot)} holds all values of \mathcmd{x} for which the upper bound function \upperBoundFunctionDot (Definition \ref{def:upper_bound_function}) is larger than the maximum observed value \mathcmd{\mathcmd{y}^{+}} and is defined by
    \begin{align*}
        \mathcmd{\upperBoundRegion = \{ \mathcmd{x} \in \domainOperator(\mathcmd{f}) ~:~ \upperBoundFunction > \observationMaxSample \} }.
    \end{align*}
\end{defn}
Note that by definition \mathcmd{\mathcal{U}} does not include any observations in \mathcmd{\mathcal{O}} since $\upperBoundFunction = \mathcmd{\mathcmd{y}_{\mathcmd{s}}} \leq \mathcmd{\mathcmd{y}^{+}}$ for all \mathcmd{\mathcmd{x}_{\mathcmd{s}}}.

\begin{lemma}
    \label{lemma:optimumInU}
    The upper bound region \mathcmd{\mathcal{U}} includes the optimum sample \mathcmd{\mathcmd{x}^*} if it has not been observed.
    Formally, if $\mathcmd{\mathcmd{x}^{+}} \neq \mathcmd{\mathcmd{x}^*}$, then $\mathcmd{\mathcmd{x}^*} \in \mathcmd{\mathcal{U}}$.
    \begin{proof}
        By Definition \ref{def:upper_bound_function} and \cref{eq:LipschitzConstant}, the value of the upper bound function at the optimal sample is larger than the optimal value.
        Formally, $\bar{\mathcmd{f}}(\mathcmd{\mathcmd{x}^*} | \mathcmd{\mathcal{O}}) \geq \mathcmd{\mathcmd{y}^*}$.
        If the optimal sample has not been observed ($\mathcmd{\mathcmd{x}^*} \neq \mathcmd{\mathcmd{x}^{+}}$), then $\mathcmd{\mathcmd{y}^*} > \mathcmd{\mathcmd{y}^{+}}$.
        Consequently, $\bar{\mathcmd{f}}(\mathcmd{\mathcmd{x}^*} | \mathcmd{\mathcal{O}}) \geq \mathcmd{\mathcmd{y}^{+}}$.
        Therefore, by definition of the upper bound region (Definition \ref{def:upperBoundRegion}) the optimal sample is included ($\mathcmd{\mathcmd{x}^*} \in \mathcmd{\mathcal{U}}$).
    \end{proof}
\end{lemma}

\structure{how can we relate $U$ and $S$?}
Based on the definition of the upper bound region, the optimum inclusion is satisfied when the set inclusion $\mathcmd{\mathcal{U}} \subseteq \mathcmd{\mathcal{S}}$ holds.
The conditions for the set inclusion are given in two parts.
Firstly, in Lemma \ref{lemma:x_in_searchRegion} it is shown that if $\meanFunction + \mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} \geq \mathcmd{\mathcmd{y}^{+}} + \mathcmd{\xi}$ and $\mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} > 0$ then $\mathcmd{x} \in \mathcmd{\mathcal{S}}$.
Secondly, by Definition \ref{def:upperBoundRegion}, if $\upperBoundFunction > \mathcmd{\mathcmd{y}^{+}}$ then $\mathcmd{x} \in \mathcmd{\mathcal{U}}$.
By combining both conditions, we find $\mathcmd{\mathcal{U}} \subseteq \mathcmd{\mathcal{S}}$ if $\meanFunction + \mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} \geq \upperBoundFunction + \mathcmd{\xi}$ for all $\mathcmd{x} \in \domainOperator(\mathcmd{f})$, as shown in Lemma \ref{lemma:setInclusion}.

\newcommand{\w}{
    \mathcmd{
        w(\mathcmd{x}, \mathcmd{\xi} | \mathcmd{\mathcal{O}})
    }
}
\newcommand{\z}{
    \mathcmd{
        z(\mathcmd{x}, \mathcmd{\xi} | \mathcmd{\mathcal{O}})
    }
}
\newcommand{\h}{
    \mathcmd{
        h(z)
    }
}
\input{sections/Analysis/lemma_x_in_searchRegion.tex}

\begin{lemma}
    If $\meanFunction + \mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} \geq \upperBoundFunction + \mathcmd{\xi}$ for all $\mathcmd{x} \in \domainOperator(\mathcmd{f})$, then $\mathcmd{\mathcal{U}} \subseteq \mathcmd{\mathcal{S}}$.
    \begin{proof}
        As shown in Lemma \ref{lemma:x_in_searchRegion}, if $\meanFunction + \mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} \geq \mathcmd{\mathcmd{y}^{+}} + \mathcmd{\xi}$ and $\mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} > 0$, then $\mathcmd{x} \in \mathcmd{\mathcal{S}}$.
        By definition of \mathcmd{\mathcal{U}}, for all $\mathcmd{x} \in \mathcmd{\mathcal{U}}$ we find $\upperBoundFunction > \mathcmd{\mathcmd{y}^{+}}$.
        Additionally, for all $\mathcmd{x} \in \mathcmd{\mathcal{U}}$ we find $\mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} > 0$, since $\mathcmd{\varianceFunctionBase(\mathcmd{x} | \observationSet)} = 0$ only if $\mathcmd{x} = \observationFunctionInput$ and $\observationFunctionInput \notin \mathcmd{\mathcal{U}}$.
        Therefore, if $\meanFunction + \mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} \geq \upperBoundFunction + \mathcmd{\xi}$, then $\mathcmd{x} \in \mathcmd{\mathcal{S}}$ for all $\mathcmd{x} \in \mathcmd{\mathcal{U}}$.
        In conclusion, $\mathcmd{\mathcal{U}} \subseteq \mathcmd{\mathcal{S}}$.
    \end{proof}
\label{lemma:setInclusion}
\end{lemma}

\structure{how to ensure the inequality holds?}
Note that \meanFunctionDot and \mathcmd{\deviationFunctionBase(\cdot | \observationSet)} depend on the kernel and \upperBoundFunctionDot depends on the Lipschitz constant.
This raises the question of which (type of) kernel is capable of explicitly associating its mean function and variance function to the Lipschitz constant.

\structure{Use Markovian kernel for its simplicity/similarity to upper bound function}
An answer can be found in the work of \citeA{Ding2018}.
\citeA{Ding2018} introduced the Markovian class kernels as kernels that possess a Markovian posterior distribution.
A Markovian posterior distribution for a certain input only depends on the observations surrounding that input.
This property is beneficial as the upper bound function, which is directly related to the Lipschitz constant, possesses the same property.
An additional benefit of this class of kernels is that the elements of $\mathcmd{\bm{K}}^{-1}(\mathcmd{\mathcal{O}})$ can be expressed analytically.
This removes the need of inversion of a matrix of which the size grows with the number of observations, since $\mathcmd{\bm{K}}(\mathcmd{\mathcal{O}}) \in \mathcmd{\mathds{R}}^{\mathcmd{S} \times \mathcmd{S}}$.
As noted by \citeA{Rasmussen2004}, this inversion is considered a major restriction to the practical application of \BO.
\noindent
In general, a Markovian class kernel is defined by
\begin{align*}
    \mathcmd{\kappa}(\mathcmd{\functionInputBase_\mathcmd{i}}, \mathcmd{\functionInputBase_\mathcmd{j}}) = \mathcmd{\lambda}^2 \Big( p(\mathcmd{\functionInputBase_\mathcmd{i}})g(\mathcmd{\functionInputBase_\mathcmd{j}}) \mathds{I}_{\mathcmd{\functionInputBase_\mathcmd{i}} \leq \mathcmd{\functionInputBase_\mathcmd{j}}} +p(\mathcmd{\functionInputBase_\mathcmd{j}})g(\mathcmd{\functionInputBase_\mathcmd{i}}) \mathds{I}_{\mathcmd{\functionInputBase_\mathcmd{i}} > \mathcmd{\functionInputBase_\mathcmd{j}}} \Big)
\end{align*}
for some function $p(\cdot)$ and $g(\cdot)$, where $\mathds{I}(\cdot)$ is the indicator function and \mathcmd{\lambda} is the kernel scale parameter.
The mean function \meanFunctionObservationDot and the variance function $\mathcmd{\deviationFunctionBase^2}_\mathcmd{s}(\cdot | \mathcmd{\mathcal{O}})$ of the posterior on the interval between observations for a kernel of this class is defined by,
\begin{align}
    \meanFunctionObservation &=
        \mathcmd{\bm{\kappa}}^\mathrm{T}_\mathcmd{s} (\mathcmd{x}, \mathcmd{\mathcal{O}})
        \mathcmd{\bm{K}}^{-1}_\mathcmd{s} (\mathcmd{\mathcal{O}})
        \mathcmd{\bm{y}}_\mathcmd{s} (\mathcmd{\mathcal{O}})
    \label{eq:mean_function_interval}
    \\
    \mathcmd{\deviationFunctionBase^2}_\mathcmd{s}(\mathcmd{x} | \mathcmd{\mathcal{O}}) &=
    \mathcmd{\kappa}\left( \mathcmd{x}, \mathcmd{x} \right) - \mathcmd{\bm{\kappa}}_\mathcmd{s}^\mathrm{T} (\mathcmd{x}, \mathcmd{\mathcal{O}})
    \mathcmd{\bm{K}}^{-1}_\mathcmd{s}(\mathcmd{\mathcal{O}})
    \mathcmd{\bm{\kappa}}_\mathcmd{s} (\mathcmd{x}, \mathcmd{\mathcal{O}})
    \label{eq:variance_function_interval}
\end{align}
for $\mathcmd{x} \in [\functionInputObservationLeft,  \functionInputObservation]$, where
\begin{align*}
    \mathcmd{\bm{\kappa}}_\mathcmd{s}(\mathcmd{x}, \mathcmd{\mathcal{O}})&=
    \begin{bmatrix}
        \mathcmd{\kappa}\bracket{\mathcmd{x}_1, \mathcmd{x}} &
        \dots &
        \mathcmd{\kappa}\bracket{\functionInputObservationLeft, \mathcmd{x}} &
        \mathcmd{\kappa}\bracket{\functionInputObservation, \mathcmd{x}} &
        \mathcmd{\kappa}\bracket{\functionInputObservationRight, \mathcmd{x}} &
        \dots &
        \mathcmd{\kappa}\bracket{\mathcmd{x}_\mathcmd{S}, \mathcmd{x}}
    \end{bmatrix}^\mathrm{T},
    \\
    \mathcmd{\bm{y}}_\mathcmd{s}(\mathcmd{\mathcal{O}}) &=
    \begin{bmatrix}
        \mathcmd{y}_1 &
        \dots &
        \functionOutputObservationLeft &
        \functionOutputObservation &
        \functionOutputObservationRight &
        \dots &
        \mathcmd{y}_\mathcmd{S}
    \end{bmatrix}^\mathrm{T},
\end{align*}
and $\mathcmd{\bm{K}}^{-1}_\mathcmd{s}(\mathcmd{\mathcal{O}})$ is a tridiagonal matrix of appropriate dimensions where the tridiagonal elements of $\mathcmd{\bm{K}}_\mathcmd{s}^{-1}(\mathcmd{\mathcal{O}})$, for $\mathcmd{S} \geq 3$ and if $\mathcmd{\bm{K}}_\mathcmd{s}(\mathcmd{\mathcal{O}})$ is nonsingular, are given by
\begin{align*}
    (\mathcmd{\bm{K}}_\mathcmd{s}^{-1}(\mathcmd{\mathcal{O}}))_{\mathcmd{s},\mathcmd{s}}
    =
    \begin{cases}
        \frac{
            \mathcmd{\lambda}^{-2}
            p(\mathcmd{x}_2)
        }{
            p(\mathcmd{x}_1)
            \bracketbig{
                p(\mathcmd{x}_2) g(\mathcmd{x}_1) - p(\mathcmd{x}_1) g(\mathcmd{x}_2)
            }
        }, &\text{if } \mathcmd{s} = 1,
        \\[12pt]
        \frac{
            \mathcmd{\lambda}^{-2}
            \bracketbig{
                p(\mathcmd{x}_{\mathcmd{s}+1}) g(\mathcmd{x}_{\mathcmd{s}-1}) - p(\mathcmd{x}_{\mathcmd{s}-1})g(\mathcmd{x}_{\mathcmd{s}+1})
            }
        }{
            \bracketbig{
                p(\mathcmd{x}_{\mathcmd{s}}) g(\mathcmd{x}_{\mathcmd{s}-1}) - p(\mathcmd{x}_{\mathcmd{s}-1}) g(\mathcmd{x}_{\mathcmd{s}})
            }
            \bracketbig{
                p(\mathcmd{x}_{\mathcmd{s}+1})  g(\mathcmd{x}_{\mathcmd{s}}) - p(\mathcmd{x}_{\mathcmd{s}}) g(\mathcmd{x}_{\mathcmd{s}+1})
            }
        }, &\text{if } \mathcmd{s} \in \{ 2, \dots, \mathcmd{S} - 1\},
        \\[12pt]
        \frac{
            \mathcmd{\lambda}^{-2}
            g(\mathcmd{x}_{\mathcmd{S}-1})
        }{
            g(\mathcmd{x}_{\mathcmd{S}})
            \bracketbig{
                p(\mathcmd{x}_{\mathcmd{S}}) g(\mathcmd{x}_{\mathcmd{S}-1}) - p(\mathcmd{x}_{\mathcmd{S} - 1}) g(\mathcmd{x}_{\mathcmd{S}})
            }
        }, &\text{if } \mathcmd{s} = \mathcmd{S},
    \end{cases}
\end{align*}
and
\begin{align*}
    (\mathcmd{\bm{K}}_\mathcmd{s}^{-1}(\mathcmd{\mathcal{O}}))_{\mathcmd{s}-1,\mathcmd{s}}
    =
    (\mathcmd{\bm{K}}_\mathcmd{s}^{-1}(\mathcmd{\mathcal{O}}))_{\mathcmd{s},\mathcmd{s}-1}
    =
    \frac{
        -\mathcmd{\lambda}^{-2}
    }{
        \bracketbig{
            p(\mathcmd{x}_{\mathcmd{s}}) g(\mathcmd{x}_{\mathcmd{s}-1}) - p(\mathcmd{x}_{\mathcmd{s}-1}) g(\mathcmd{x}_{\mathcmd{s}})
        }
    }
    ,~\mathcmd{s} = 2,\dots,\mathcmd{S}.
\end{align*}
All other elements of $\mathcmd{\bm{K}}_\mathcmd{s}^{-1}(\mathcmd{\mathcal{O}})$ are equal to zero.

\structure{which kernel and what parameters to choose?}
Next, we show that for the Dirichlet kernel, as introduced by \citeA{Ding2018}, the inequality of Lemma \ref{lemma:setInclusion} will hold for all observations if the kernel scale is chosen appropriately.
This kernel is selected over other Markovian class kernels because of to its simplicity.
The Dirichlet kernel defined by
\begin{align}
    \label{eq:dirichlet}
    \dirichletKernelLong
\end{align}
for $\mathcmd{\functionInputBase_\mathcmd{i}}, \mathcmd{\functionInputBase_\mathcmd{j}} \in [0, 1]$.
Note that for \mathcmd{\kernel_{\mathrm{d}}}, we find that $p(\mathcmd{x}) = \mathcmd{x}$ and $g(\mathcmd{x}) = (1 - \mathcmd{x})$.
The mean function, given by \cref{eq:mean_function_interval}, and the variance function, given by \cref{eq:variance_function_interval}, corresponding to the Dirichlet kernel in the interval $[\functionInputObservationLeft, \functionInputObservation]$ can be written as
\begin{align}
    \meanFunctionObservation &=
    \frac{\mathcmd{y}_{\mathcmd{s}-1}(\functionInputObservation - x) + \mathcmd{y}_{\mathcmd{s}}(x - \functionInputObservationLeft)}{\functionInputObservation - \functionInputObservationLeft},
    \label{eq:mean_function_interval_Dirichlet}
\end{align}
\begin{align}
    \mathcmd{\deviationFunctionBase^2}_\mathcmd{s}(\mathcmd{x} | \mathcmd{\mathcal{O}})
    &=
        \mathcmd{\lambda}^2
        \frac{
            -(\functionInputObservation - \mathcmd{x})
            (\functionInputObservationLeft - \mathcmd{x})
        }{
            \functionInputObservation - \functionInputObservationLeft
        }.
    \label{eq:variance_function_interval_Dirichlet}
\end{align}
The derivation of \cref{eq:mean_function_interval_Dirichlet,eq:variance_function_interval_Dirichlet} can be found in \cref{section:dirichlet_interval_functions}.
Note that both the mean function and the variance function on the interval $[\functionInputObservationLeft, \functionInputObservation]$ only depend on the observations ($\mathcmd{O}_{\mathcmd{s}-1}$, $\mathcmd{O}_{\mathcmd{s}}$) at the boundaries of the interval.

\structure{Sufficient to proof for a single segment}
Based on \cref{eq:mean_function_interval_Dirichlet,eq:variance_function_interval_Dirichlet}, the inequality of Lemma \ref{lemma:setInclusion} holds if
\begin{align}
    \meanFunctionObservation  + \deviationFunctionObservation
    &\geq
    \upperBoundFunction + \mathcmd{\xi}
    \text{ for }\mathcmd{x} \in [\functionInputObservationLeft, \functionInputObservation]
    \label{eq:inequality_observation}
\end{align}
for all $\mathcmd{s} \in \{1, \dots, \mathcmd{S}\}$, given $\mathcmd{x}_1 = 0$ and $\mathcmd{x}_\mathcmd{S} = 1$.
In other words, by using the Dirichlet kernel, instead of analyzing the inequality in Lemma \ref{lemma:setInclusion} over the entire domain of the function, it is sufficient to analyze \cref{eq:inequality_observation} on the intervals between the observations.

\structure{Theorem (to conclude the section)}
Now that the acquisition function and the kernel have been selected, we need to find their parameters (\mathcmd{\xi} and \mathcmd{\lambda}) such that the inequality of \cref{eq:inequality_observation} holds.
These parameters can be appropriately chosen based on the Lipschitz constant as follows:

\input{sections/Analysis/theorem_kernel_scale.tex}

According to Theorem \ref{theorem:kernel_scale}, if $\mathcmd{\lambda} \geq \lipschitzConstantIndexed$ we find $\meanFunction + \mathcmd{\deviationFunctionBase(\mathcmd{x} | \observationSet)} \geq \upperBoundFunction$ for all \mbox{$\mathcmd{x} \in \domainOperator(\mathcmd{f})$}.
Applying Lemma \ref{lemma:setInclusion} and setting $\mathcmd{\xi} = 0$, yields $\mathcmd{\mathcal{U}} \subseteq \mathcmd{\mathcal{S}}$ where $\mathcmd{\mathcmd{x}^*} \in \mathcmd{\mathcal{S}}$ for all observations.
Subsequently, the Bayesian optimization will converge to the global optimum of the function.
Note that for all functions without normalized domain, the Lipschitz constant should be scaled according to the scaling required for normalization of the domain.

\section{Distributed Bayesian Algorithm}
\label{section:distributed_bayesian_algorithm}

\newcommand{\algrule}[1][.2pt]{\par\vskip.5\baselineskip\hrule height #1\par\vskip.5\baselineskip}

\begin{algorithm}

\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{$\mathcmd{\Parent_{\mathcmd{i}}}, \mathcmd{\PseudoParent_{\mathcmd{i}}}, \mathcmd{\Child_{\mathcmd{i}}}, \mathcmd{\PseudoChild_{\mathcmd{i}}}, \functionSetLocalIndexed, \functionSetSharedIndexed, \localVariablesIndexed, \mathcmd{\kappa}$}
\Output{$\hat{\mathcmd{\sigma}}_{\localVariablesIndexed}$}

\SetKwFor{ForEach}{foreach}{do}{}
\SetKwFor{When}{when received}{}{}
\SetKwBlock{Initialization}{Initialization}{}
\SetKwData{Root}{root agent}
\SetKwData{Leaf}{leaf agent}
\SetKwFunction{FnSample}{computeOptimalSample(\mathcmd{\kappa})}
\SetKwFunction{FnStoreUtility}{storeUtility(\optimumMessageToParent, \sampleMessageIndexed)}

\SetKwFunction{FnRetrieveBestSample}{retrieveOptimalLocalSample(\mathcmd{\hat{\sampleMessage}_\agentIndexParent})}
\SetKwFunction{FnRetrieveBestUtility}{retrieveOptimalUtility()}

\SetKwFunction{FnSendSample}{send(\mathcmd{\agent_\agentIndexChild}, \sampleMessageIndexed)}
\SetKwFunction{FnSendSampleInitial}{send(\mathcmd{\agent_\agentIndexChild}, \sampleMessageIndexed)}

\SetKwFunction{FnSendSampleBest}{send(\mathcmd{\agent_\agentIndexChild}, \mathcmd{\hat{\sampleMessage}_\mathcmd{i}})}
\SetKwFunction{FnSendSampleBestInitial}{send(\mathcmd{\agent_\agentIndexChild}, \mathcmd{\hat{\sampleMessage}_\mathcmd{i}})}

\SetKwFunction{FnSendUtility}{send(\mathcmd{\Parent_{\mathcmd{i}}}, \optimumMessageToParent)}

\SetKwFunction{SubOptimizeLocal}{optimizeLocalVariables(\mathcmd{\sampleMessage_\agentIndexParent})}
\SetKwFunction{SubOptimizeLocalRoot}{optimizeLocalVariables($\emptyset$)}

\SetKwFunction{SubSamplePhase}{createSampleMessage(\mathcmd{\sampleMessage_\agentIndexParent})}
\SetKwFunction{SubSamplePhaseRoot}{createSampleMessage($\emptyset$)}

\SetKwFunction{SubUtilityPhase}{calculateUtility(\optimumMessageBestChild)}
\SetKwFunction{SubUtilityPhaseLeaf}{calculateUtility()}

\SetKwFunction{SubFinalPhase}{processFinal(\mathcmd{\hat{\sampleMessage}_\agentIndexParent})}
\SetKwFunction{SubFinalPhaseRoot}{processFinal($\emptyset$)}

\SetKwFunction{SubCalculateUtility}{calculateUtility(\sampleMessageIndexed)}
\SetKwFunction{SubGetChildUtility}{getChildUtility(\sampleMessageIndexed)}

\Initialization{
    \uIf{\Root}{
    \While{not threshold reached}{
        $\optimumMessageToParent \coloneqq \SubOptimizeLocalRoot$\;
    }
        \SubFinalPhaseRoot\;
    }
}
\When{\mathcmd{\mathrm{\bf sample}} \mathcmd{\sampleMessage_\agentIndexParent} from parent \mathcmd{\Parent_{\mathcmd{i}}}}{

    \While{not threshold reached}{
        $\optimumMessageToParent \coloneqq \SubOptimizeLocal$\;
    }
    \FnSendUtility\;
}

\When{\mathcmd{\mathrm{\bf final}} \mathcmd{\hat{\sampleMessage}_\agentIndexParent} from parent \mathcmd{\Parent_{\mathcmd{i}}}}{
    \SubFinalPhase\;
}

\caption{\DBayLong for agent \mathcmd{\agent_\agentIndex}}
\label{algorithm:DBAY}
\end{algorithm}

\begin{function}
    \SetKwFor{ForEach}{foreach}{do}{}
    \SetKwFor{When}{when received}{}{}
    \SetKwBlock{Initialization}{Initialization}{}
    \SetKwData{Root}{root agent}
    \SetKwData{Leaf}{leaf agent}
    \SetKwFunction{FnSample}{computeOptimalSample(\mathcmd{\kappa})}
    \SetKwFunction{FnStoreUtility}{storeUtility(\optimumMessageToParent, \sampleMessageIndexed)}

    \SetKwFunction{FnRetrieveBestSample}{retrieveOptimalLocalSample(\mathcmd{\hat{\sampleMessage}_\agentIndexParent})}
    \SetKwFunction{FnRetrieveBestUtility}{retrieveOptimalUtility()}

    \SetKwFunction{FnSendSample}{send(\mathcmd{\agent_\agentIndexChild}, \sampleMessageIndexed)}
    \SetKwFunction{FnSendSampleInitial}{send(\mathcmd{\agent_\agentIndexChild}, \sampleMessageIndexed)}

    \SetKwFunction{FnSendSampleBest}{send(\mathcmd{\agent_\agentIndexChild}, \mathcmd{\hat{\sampleMessage}_\mathcmd{i}})}
    \SetKwFunction{FnSendSampleBestInitial}{send(\mathcmd{\agent_\agentIndexChild}, \mathcmd{\hat{\sampleMessage}_\mathcmd{i}})}

    \SetKwFunction{FnSendUtility}{send(\mathcmd{\Parent_{\mathcmd{i}}}, \optimumMessageToParent)}
    \SetKwFunction{SubOptimizeLocal}{optimizeLocalVariables(\mathcmd{\sampleMessage_\agentIndexParent})}
    \SetKwFunction{SubOptimizeLocalRoot}{optimizeLocalVariables($\emptyset$)}

    \SetKwFunction{SubSamplePhase}{createSampleMessage(\mathcmd{\sampleMessage_\agentIndexParent})}
    \SetKwFunction{SubSamplePhaseRoot}{createSampleMessage($\emptyset$)}

    \SetKwFunction{SubUtilityPhase}{calculateUtility(\optimumMessageBestChild)}
    \SetKwFunction{SubUtilityPhaseLeaf}{calculateUtility()}

    \SetKwFunction{SubFinalPhase}{processFinal(\mathcmd{\hat{\sampleMessage}_\agentIndexParent})}
    \SetKwFunction{SubFinalPhaseRoot}{processFinal($\emptyset$)}

    \SetKwFunction{SubCalculateUtility}{calculateUtility(\sampleMessageIndexed)}
    \SetKwFunction{SubGetChildUtility}{getChildUtility(\sampleMessageIndexed)}

    \SetKwProg{ProgSubOptimizeLocal}{Function}{}{}
    \ProgSubOptimizeLocal{\SubOptimizeLocal}{
    $\mathcmd{\sigma}_{\localVariablesIndexed} \coloneqq \FnSample$\;
    \sampleMessageIndexedOtherLongInAlgorithm\;
    $\optimumMessageToParent \coloneqq \SubCalculateUtility$\;
    \KwRet \optimumMessageToParent\;
    }
\end{function}

\clearpage
\begin{function}
    \SetKwFor{ForEach}{foreach}{do}{}
    \SetKwFor{When}{when received}{}{}
    \SetKwBlock{Initialization}{Initialization}{}
    \SetKwData{Root}{root agent}
    \SetKwData{Leaf}{leaf agent}
    \SetKwFunction{FnSample}{computeOptimalSample(\mathcmd{\kappa})}
    \SetKwFunction{FnStoreUtility}{storeUtility(\optimumMessageToParent, \sampleMessageIndexed)}

    \SetKwFunction{FnRetrieveBestSample}{retrieveOptimalLocalSample(\mathcmd{\hat{\sampleMessage}_\agentIndexParent})}
    \SetKwFunction{FnRetrieveBestUtility}{retrieveOptimalUtility()}

    \SetKwFunction{FnSendSample}{send(\mathcmd{\agent_\agentIndexChild}, \sampleMessageIndexed)}
    \SetKwFunction{FnSendSampleInitial}{send(\mathcmd{\agent_\agentIndexChild}, \sampleMessageIndexed)}

    \SetKwFunction{FnSendSampleBest}{send(\mathcmd{\agent_\agentIndexChild}, \mathcmd{\hat{\sampleMessage}_\mathcmd{i}})}
    \SetKwFunction{FnSendSampleBestInitial}{send(\mathcmd{\agent_\agentIndexChild}, \mathcmd{\hat{\sampleMessage}_\mathcmd{i}})}

    \SetKwFunction{FnSendUtility}{send(\mathcmd{\Parent_{\mathcmd{i}}}, \optimumMessageToParent)}
    \SetKwFunction{SubOptimizeLocal}{optimizeLocalVariables(\mathcmd{\sampleMessage_\agentIndexParent})}
    \SetKwFunction{SubOptimizeLocalRoot}{optimizeLocalVariables($\emptyset$)}

    \SetKwFunction{SubSamplePhase}{createSampleMessage(\mathcmd{\sampleMessage_\agentIndexParent})}
    \SetKwFunction{SubSamplePhaseRoot}{createSampleMessage($\emptyset$)}

    \SetKwFunction{SubUtilityPhase}{calculateUtility(\optimumMessageBestChild)}
    \SetKwFunction{SubUtilityPhaseLeaf}{calculateUtility()}

    \SetKwFunction{SubFinalPhase}{processFinal(\mathcmd{\hat{\sampleMessage}_\agentIndexParent})}
    \SetKwFunction{SubFinalPhaseRoot}{processFinal($\emptyset$)}

    \SetKwFunction{SubCalculateUtility}{calculateUtility(\sampleMessageIndexed)}
    \SetKwFunction{SubGetChildUtility}{getChildUtility(\sampleMessageIndexed)}

    \SetKwProg{ProgSubCalculateUtility}{Function}{}{}
    \ProgSubCalculateUtility{\SubCalculateUtility}{
    \optimumMessageLocalLongInAlgorithm\;
    \uIf{$\mathcmd{\Child_{\mathcmd{i}}} \neq \emptyset$}{
        $\optimumMessageBestChild \coloneqq \SubGetChildUtility$\;
        \optimumMessageToParentLongInAlgorithm\;
    } \uElse{
        $\optimumMessageToParent \coloneqq \optimumMessageLocal$\;
    }
    \FnStoreUtility\;
        \KwRet \optimumMessageToParent\;
    }
\end{function}

\begin{function}
    \SetKwFor{ForEach}{foreach}{do}{}
    \SetKwFor{When}{when received}{}{}
    \SetKwBlock{Initialization}{Initialization}{}
    \SetKwData{Root}{root agent}
    \SetKwData{Leaf}{leaf agent}
    \SetKwFunction{FnSample}{computeOptimalSample(\mathcmd{\kappa})}
    \SetKwFunction{FnStoreUtility}{storeUtility(\optimumMessageToParent, \sampleMessageIndexed)}

    \SetKwFunction{FnRetrieveBestSample}{retrieveOptimalLocalSample(\mathcmd{\hat{\sampleMessage}_\agentIndexParent})}
    \SetKwFunction{FnRetrieveBestUtility}{retrieveOptimalUtility()}

    \SetKwFunction{FnSendSample}{send(\mathcmd{\agent_\agentIndexChild}, \sampleMessageIndexed)}
    \SetKwFunction{FnSendSampleInitial}{send(\mathcmd{\agent_\agentIndexChild}, \sampleMessageIndexed)}

    \SetKwFunction{FnSendSampleBest}{send(\mathcmd{\agent_\agentIndexChild}, \mathcmd{\hat{\sampleMessage}_\mathcmd{i}})}
    \SetKwFunction{FnSendSampleBestInitial}{send(\mathcmd{\agent_\agentIndexChild}, \mathcmd{\hat{\sampleMessage}_\mathcmd{i}})}

    \SetKwFunction{FnSendUtility}{send(\mathcmd{\Parent_{\mathcmd{i}}}, \optimumMessageToParent)}
    \SetKwFunction{SubOptimizeLocal}{optimizeLocalVariables(\mathcmd{\sampleMessage_\agentIndexParent})}
    \SetKwFunction{SubOptimizeLocalRoot}{optimizeLocalVariables($\emptyset$)}

    \SetKwFunction{SubSamplePhase}{createSampleMessage(\mathcmd{\sampleMessage_\agentIndexParent})}
    \SetKwFunction{SubSamplePhaseRoot}{createSampleMessage($\emptyset$)}

    \SetKwFunction{SubUtilityPhase}{calculateUtility(\optimumMessageBestChild)}
    \SetKwFunction{SubUtilityPhaseLeaf}{calculateUtility()}

    \SetKwFunction{SubFinalPhase}{processFinal(\mathcmd{\hat{\sampleMessage}_\agentIndexParent})}
    \SetKwFunction{SubFinalPhaseRoot}{processFinal($\emptyset$)}

    \SetKwFunction{SubCalculateUtility}{calculateUtility(\sampleMessageIndexed)}
    \SetKwFunction{SubGetChildUtility}{getChildUtility(\sampleMessageIndexed)}

    \SetKwProg{ProgSubGetChildUtility}{Function}{}{}
    \ProgSubGetChildUtility{\SubGetChildUtility}{
    \ForEachSingleLine{$\mathcmd{\agent_\agentIndexChild} \in \mathcmd{\Child_{\mathcmd{i}}}$}{
        \FnSendSample\;
    }
    \When{\optimumMessageFromChild from all $\mathcmd{\agent_\agentIndexChild} \in \mathcmd{\Child_{\mathcmd{i}}}$}{
        \optimumMessageIndexedChildrenInAlgorithm\;
    }
        \KwRet \optimumMessageBestChild\;
    }
\end{function}

\begin{function}
    \SetKwFor{ForEach}{foreach}{do}{}
    \SetKwFor{When}{when received}{}{}
    \SetKwBlock{Initialization}{Initialization}{}
    \SetKwData{Root}{root agent}
    \SetKwData{Leaf}{leaf agent}
    \SetKwFunction{FnSample}{computeOptimalSample(\mathcmd{\kappa})}
    \SetKwFunction{FnStoreUtility}{storeUtility(\optimumMessageToParent, \sampleMessageIndexed)}

    \SetKwFunction{FnRetrieveBestSample}{retrieveOptimalLocalSample(\mathcmd{\hat{\sampleMessage}_\agentIndexParent})}
    \SetKwFunction{FnRetrieveBestUtility}{retrieveOptimalUtility()}

    \SetKwFunction{FnSendSample}{send(\mathcmd{\agent_\agentIndexChild}, \sampleMessageIndexed)}
    \SetKwFunction{FnSendSampleInitial}{send(\mathcmd{\agent_\agentIndexChild}, \sampleMessageIndexed)}

    \SetKwFunction{FnSendSampleBest}{send(\mathcmd{\agent_\agentIndexChild}, \mathcmd{\hat{\sampleMessage}_\mathcmd{i}})}
    \SetKwFunction{FnSendSampleBestInitial}{send(\mathcmd{\agent_\agentIndexChild}, \mathcmd{\hat{\sampleMessage}_\mathcmd{i}})}

    \SetKwFunction{FnSendUtility}{send(\mathcmd{\Parent_{\mathcmd{i}}}, \optimumMessageToParent)}
    \SetKwFunction{SubOptimizeLocal}{optimizeLocalVariables(\mathcmd{\sampleMessage_\agentIndexParent})}
    \SetKwFunction{SubOptimizeLocalRoot}{optimizeLocalVariables($\emptyset$)}

    \SetKwFunction{SubSamplePhase}{createSampleMessage(\mathcmd{\sampleMessage_\agentIndexParent})}
    \SetKwFunction{SubSamplePhaseRoot}{createSampleMessage($\emptyset$)}

    \SetKwFunction{SubUtilityPhase}{calculateUtility(\optimumMessageBestChild)}
    \SetKwFunction{SubUtilityPhaseLeaf}{calculateUtility()}

    \SetKwFunction{SubFinalPhase}{processFinal(\mathcmd{\hat{\sampleMessage}_\agentIndexParent})}
    \SetKwFunction{SubFinalPhaseRoot}{processFinal($\emptyset$)}

    \SetKwFunction{SubCalculateUtility}{calculateUtility(\sampleMessageIndexed)}
    \SetKwFunction{SubGetChildUtility}{getChildUtility(\sampleMessageIndexed)}

    \SetKwProg{ProgFinalPhase}{Function}{}{}
    \ProgFinalPhase{\SubFinalPhase}{
    $\hat{\mathcmd{\sigma}}_{\localVariablesIndexed} \coloneqq \FnRetrieveBestSample$\;
    \sampleMessageBestIndexedLongInAlgorithm\;
    \ForEachSingleLine{$\mathcmd{\agent_\agentIndexChild} \in \mathcmd{\Child_{\mathcmd{i}}}$}{
        \FnSendSampleBest\;
    }
    }
\end{function}
\section{Dirichlet kernel interval functions}
\label{section:dirichlet_interval_functions}
\structure{General Markovian property}
As shown in the work of \citeA[Theorem~2]{Ding2018}, a kernel \mathcmd{\kappa}  of the Markovian class reduces the mean function \meanFunctionObservationDot and the variance function \varianceFunctionObservationDot of the posterior on the interval between observations as given in \cref{eq:mean_function_interval,eq:variance_function_interval}, respectively.
For the Dirichlet kernel as defined by \cref{eq:dirichlet}, for a normalized domain $\mathcmd{\functionInputBase_\mathcmd{i}}, \mathcmd{\functionInputBase_\mathcmd{j}} \in [0, 1]$ and the kernel scale parameter \mathcmd{\lambda}, the non-zero elements of the $\mathcmd{\bm{K}}^{-1}_\mathcmd{s}(\mathcmd{\mathcal{O}})$ matrix are given by
\begin{align*}
    (\mathcmd{\bm{K}}^{-1}_\mathcmd{s}(\mathcmd{\mathcal{O}}))_{\mathcmd{s},\mathcmd{s}}
    =
    \begin{cases}
        \mathcmd{\lambda}^{-2}
        \frac{
            \mathcmd{x}_1
        }{
            \mathcmd{x}_1
            \bracketbig{
                \mathcmd{x}_2 - \mathcmd{x}_1
            }
        }
        , &\text{if } \mathcmd{s} = 1,
        \\[12pt]
        \mathcmd{\lambda}^{-2}
        \frac{
            (\functionInputObservationRight - \functionInputObservationLeft)
        }{
            (\functionInputObservation - \functionInputObservationLeft)
            (\functionInputObservationRight - \functionInputObservation)
        }
        , &\text{if } \mathcmd{s} \in \{ 2, \dots, \mathcmd{S} - 1\},
        \\[12pt]
        \mathcmd{\lambda}^{-2}
        \frac{
            (1 - \mathcmd{x}_{\mathcmd{S}-1})
        }{
            (1 - \mathcmd{x}_\mathcmd{S})
            (\mathcmd{x}_\mathcmd{S} - \mathcmd{x}_{\mathcmd{S} - 1})
        }
        , &\text{if } \mathcmd{s} = \mathcmd{S},
    \end{cases}
\end{align*}
and
\begin{align*}
    (\mathcmd{\bm{K}}_\mathcmd{s}^{-1}(\mathcmd{\mathcal{O}}))_{\mathcmd{s}-1,\mathcmd{s}}
    =
    (\mathcmd{\bm{K}}_\mathcmd{s}^{-1}(\mathcmd{\mathcal{O}}))_{\mathcmd{s},\mathcmd{s}-1}
    =
    \frac{
        -\mathcmd{\lambda}^{-2}
    }{
        (\functionInputObservation - \functionInputObservationLeft)
    }
    ,~\mathcmd{s} = 2,\dots,\mathcmd{S}.
\end{align*}

The mean function \meanFunctionObservationDot and the variance function \varianceFunctionObservationDot for the Dirichlet kernel can be rewritten accordingly as
\begin{align}
    \meanFunctionObservation &=
    \mathcmd{\bm{\kappa}}^\mathrm{T}_\mathcmd{s} (\mathcmd{x}, \mathcmd{\mathcal{O}})
    \mathcmd{\bm{K}}^{-1}_\mathcmd{s} (\mathcmd{\mathcal{O}})
    \mathcmd{\bm{y}}_\mathcmd{s} (\mathcmd{\mathcal{O}})
    \nonumber \\
    &=
    \begin{bmatrix}
        0
        &
        \dots
        &
        0
        &
        \frac{\functionInputObservation - x}{\functionInputObservation - \functionInputObservationLeft}
        &
        \frac{x - \functionInputObservationLeft}{\functionInputObservation - \functionInputObservationLeft}
        &
        0
        &
        \dots
        &
        0
    \end{bmatrix}
    \begin{bmatrix}
        \mathcmd{y}_1 \\
        \vdots \\
        \functionOutputObservationLeftLeft \\
        \functionOutputObservationLeft \\
        \functionOutputObservation \\
        \functionOutputObservationRight \\
        \vdots \\
        \mathcmd{y}_\mathcmd{S}
    \end{bmatrix}
    \nonumber\\
    &=
    \begin{bmatrix}
        \frac{\functionInputObservation - x}{\functionInputObservation - \functionInputObservationLeft}
        &
        \frac{x - \functionInputObservationLeft}{\functionInputObservation - \functionInputObservationLeft}
    \end{bmatrix}
    \begin{bmatrix}
        \mathcmd{y}_{\mathcmd{s}-1}
        \\
        \mathcmd{y}_{\mathcmd{s}}
    \end{bmatrix}
    \nonumber\\
    &=
    \frac{\mathcmd{y}_{\mathcmd{s}-1}(\functionInputObservation - x) + \mathcmd{y}_{\mathcmd{s}}(x - \functionInputObservationLeft)}{\functionInputObservation - \functionInputObservationLeft}
\end{align}

and
\begin{align}
    \mathcmd{\deviationFunctionBase^2}_\mathcmd{s}(\mathcmd{x} | \mathcmd{\mathcal{O}}) &=
    \mathcmd{\kappa}\left( \mathcmd{x}, \mathcmd{x} \right) - \mathcmd{\bm{\kappa}}_\mathcmd{s}^\mathrm{T} (\mathcmd{x}, \mathcmd{\mathcal{O}})
    \mathcmd{\bm{K}}^{-1}_\mathcmd{s}(\mathcmd{\mathcal{O}})
    \mathcmd{\bm{\kappa}}_\mathcmd{s} (\mathcmd{x}, \mathcmd{\mathcal{O}})
    \nonumber \\
    &=
    \mathcmd{\lambda}^2 \mathcmd{x} ( 1 - \mathcmd{x} )
    -
    \begin{bmatrix}
        0
        &
        \dots
        &
        0
        &
        \frac{\functionInputObservation - x}{\functionInputObservation - \functionInputObservationLeft}
        &
        \frac{x - \functionInputObservationLeft}{\functionInputObservation - \functionInputObservationLeft}
        &
        0
        &
        \dots
        &
        0
    \end{bmatrix}
    \begin{bmatrix}
        \mathcmd{\lambda}^2 \mathcmd{x}_1 (1 - \mathcmd{x})
        \\
        \vdots
        \\
        \mathcmd{\lambda}^2 \functionInputObservationLeftLeft (1 - \mathcmd{x})
        \\
        \mathcmd{\lambda}^2 \functionInputObservationLeft (1 - \mathcmd{x})
        \\
        \mathcmd{\lambda}^2 \mathcmd{x} (1 - \functionInputObservation)
        \\
        \mathcmd{\lambda}^2 \mathcmd{x} (1 - \functionInputObservationRight)
        \\
        \vdots
        \\
        \mathcmd{\lambda}^2 \mathcmd{x} (1 - \mathcmd{x}_\mathcmd{S})
    \end{bmatrix}
    \nonumber \\
        &=
        \mathcmd{\lambda}^2
        \bracket{
            \mathcmd{x} ( 1 - \mathcmd{x} )
            -
            \bracketBig{
                \frac{\functionInputObservationLeft (1 - \mathcmd{x}) (\functionInputObservation - x)}{\functionInputObservation - \functionInputObservationLeft}
                +
                \frac{\mathcmd{x} (1 - \functionInputObservation)(x - \functionInputObservationLeft)}{\functionInputObservation - \functionInputObservationLeft}
            }
        }
    \nonumber \\
        &=
        \mathcmd{\lambda}^2
        \frac{
            -(\functionInputObservation - \mathcmd{x})
            (\functionInputObservationLeft - \mathcmd{x})
        }{
            \functionInputObservation - \functionInputObservationLeft
        }.
\end{align}

redpajama_set_name: RedPajamaArXiv