* Sonoda et al. (2022a) S. Sonoda, I. Ishikawa, and M. Ikeda.  Universality of Group Convolutional Neural Networks Based on Ridgelet Analysis on Groups.  In _Advances in Neural Information Processing Systems 35_ , pages 38680–38694, 2022a.
  * Sonoda et al. (2022b) S. Sonoda, I. Ishikawa, and M. Ikeda.  Fully-Connected Network on Noncompact Symmetric Space and Ridgelet Transform based on Helgason-Fourier Analysis.  In _Proceedings of the 39th International Conference on Machine Learning_ , volume 162, pages 20405–20422, 2022b.
  * Sonoda et al. (2023a) S. Sonoda, Y. Hashimoto, I. Ishikawa, and M. Ikeda.  Deep Ridgelet Transform: Voice with Koopman Operator Proves Universality of Formal Deep Networks.  In _Proceedings of the 2nd NeurIPS Workshop on Symmetry and Geometry in Neural Representations_ , 2023a.
  * Sonoda et al. (2023b) S. Sonoda, H. Ishi, I. Ishikawa, and M. Ikeda.  Joint Group Invariant Functions on Data-Parameter Domain Induce Universal Neural Networks.  In _Proceedings of the 2nd NeurIPS Workshop on Symmetry and Geometry in Neural Representations_ , 2023b.
  * Starck et al. (2010) J.-L. Starck, F. Murtagh, and J. M. Fadili.  The ridgelet and curvelet transforms.  In _Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity_ , pages 89–118. Cambridge University Press, 2010.
  * Suzuki (2020) T. Suzuki.  Generalization bound of globally optimal non-convex neural network training: Transportation map estimation by infinite dimensional Langevin dynamics.  In _Advances in Neural Information Processing Systems 33_ , pages 19224–19237, 2020.
  * Terras (2016) A. Terras.  _Harmonic Analysis on Symmetric Spaces—Higher Rank Spaces, Positive Definite Matrix Space and Generalizations_.  Springer New York, 2016.
  * Unser (2019) M. Unser.  A Representer Theorem for Deep Neural Networks.  _Journal of Machine Learning Research_ , 20(110):1–30, 2019.
  * Yamasaki et al. (2023) H. Yamasaki, S. Subramanian, S. Hayakawa, and S. Sonoda.  Quantum Ridgelet Transform: Winning Lottery Ticket of Neural Networks with Quantum Computation.  In _Proceedings of the 40th International Conference on Machine Learning_ , volume 202, pages 39008–39034, 2023.
  * Yarotsky (2022) D. Yarotsky.  Universal Approximations of Invariant Maps by Neural Networks.  _Constructive Approximation_ , 55:407–474, 2022.
  * Zaheer et al. (2017) M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Poczos, R. R. Salakhutdinov, and A. J. Smola.  Deep Sets.  In _Advances in Neural Information Processing Systems_ , volume 30, 2017.
  * Zhou (2020) D.-X. Zhou.  Universality of deep convolutional neural networks.  _Applied and Computational Harmonic Analysis_ , 48(2):787–794, 2020.


# Privacy-Preserving Distributed Average Consensus in Finite Time using Random
Gossip

Nicolaos E. Manitara, Apostolos I. Rikos, and Christoforos N. Hadjicostis The
work of N. E. Manitara and C. N. Hadjicostis has been supported in part by the
European Commission’s Horizon 2020 research and innovation programme under
grant agreement No 739551 (KIOS CoE). Any opinions, findings, and conclusions,
or recommendations expressed in this publication are those of the authors and
do not necessarily reflect the views of EC.Nicolaos E. Manitara and C. N.
Hadjicostis are with the Department of Electrical and Computer Engineering and
with KIOS Center of Excelence, University of Cyprus, Nicosia 1678, Cyprus.
Emails: {manitara.nicolas<EMAIL_ADDRESS>I. Rikos is with the
Division of Decision and Control Systems, KTH Royal Institute of Technology,
SE-100 44 Stockholm, Sweden. E-mail<EMAIL_ADDRESS>
###### Abstract

In this paper, we develop and analyze a gossip-based average consensus
algorithm that enables all of the components of a distributed system, each
with some initial value, to reach (approximate) average consensus on their
initial values after executing a finite number of iterations, and without
having to reveal the specific value they contribute to the average
calculation. We consider a fully-connected (undirected) network in which each
pair of components (nodes) can be randomly selected to perform pairwise
standard gossip averaging of their values, and propose an enhancement that can
be followed by each node that does not want to reveal its initial value to
other (curious) nodes. We assume that curious nodes try to identify the
initial values of other nodes but do not interfere in the computation in any
other way; however, as a worst-case assumption, curious nodes are allowed to
collaborate arbitrarily and are assumed to know the privacy-preserving
strategy (but not the actual parameters chosen by the nodes that want to
preserve their privacy). We characterize precisely conditions on the
information exchange that guarantee privacy-preservation for a specific node.
The protocol also provides a criterion that allows the nodes to determine, in
a distributed manner (while running the enhanced gossip protocol), when to
terminate their operation because approximate average consensus has been
reached, i.e., all nodes have obtained values that are within a small distance
from the exact average of their initial values.

##  I Introduction

In distributed systems and networks, it is often necessary for all or some of
the components (nodes) to calculate a function of certain parameters that we
refer to as initial values. When all nodes calculate the average of their
initial values, they are said to reach average consensus. Average consensus
(and more generally consensus) has received a lot of attention from the
control community due to its usage in various emerging applications, including
wireless smart meters (where all nodes have to agree on the average power
demand or consumption of the network [1]), and multi-agent systems (where all
agents communicate with each other in order to coordinate their direction or
speed [2]). Over the last few decades, a variety of distributed algorithms for
calculating different functions of these initial values have been proposed by
the control, communication, and computer science communities [3, 4, 5, 6].

One approach to consensus is based on a linear iterative strategy, where each
node in the network repeatedly updates its value to be a weighted sum of its
own previous value and the values of its neighbors. If the network topology
satisfies certain conditions, the weights for the linear iteration can be
chosen such that all of the nodes in the network converge (asymptotically) to
the same function of the initial values, which under additional conditions on
the weights, can be the average [7, 8]. Another popular approach to consensus
is a gossip-based strategy [9] which allows the nodes to calculate the average
value of the network without having to follow a predefined routing of the
information that needs to be transmitted/received to/from other nodes in the
network. At each iteration, information is exchanged between a randomly
selected pair of neighboring nodes, and then this information is processed by
both nodes to compute a local pairwise average.

Both of the methods described above are asymptotic and typically do not
consider privacy issues that may arise while the network is reaching
consensus. This paper addresses the topic of privacy-preserving average
consensus in finite time, which has received some recent attention.
Specifically, an anonymization transform using random offsets on the initial
values was proposed for a cooperative111In a cooperative network, all the
nodes follow the predefined strategy, without deviating in any way [10].
wireless network in [10]. In this approach, each node that would like to
protect its privacy (i.e., does not want to reveal its initial value) adds a
random offset value to its initial value, ensuring in this way that its true
initial value will not be exposed through the values exchanged in the network.
Similarly, in the privacy-preserving strategies proposed in [11] nodes inject
uncorrelated noise into the exchanged messages so that the data associated to
a particular node cannot be inferred by a curious node during the execution of
the algorithm. However, note that the average value is not obtained exactly in
[11] leading to notions of differential privacy [12]. The injection of
correlated noise at each time step and for a finite period of time was
proposed in [13] and guarantees convergence to the exact average. In [14], the
nodes asymptotically subtract the initial offset values they added in the
computation, while in [15] each node masks its initial state with an offset
such that the sum of the offsets of each node is zero, thus guaranteeing
convergence to the average. In [16, 17], the nodes inject an initial quantized
offset to their initial values. During the operation of the algorithms, the
nodes are able to subtract the initial offset values they added, which
guarantees convergence in finite time to the exact average. In [18] the nodes
achieve privacy-preserving distributed average consensus using the principle
of additive secret sharing where each node replaces its initial value with
another obfuscated value by subtracting and adding random numbers. Finally,
another approach that guarantees privacy-preservation is via homomorphic
encryption [19, 20, 21]. However, this approach requires the existence of
trusted nodes and imposes heavier computational requirements on the nodes.

Our own previous work in [13] describes a privacy-preserving protocol that is
a variation of the standard linear iteration in [7] that is used in the
absence of privacy requirements under undirected communication topologies
using weights that form a doubly stochastic matrix. The main enhancement is
that, at each time-step, each node following the protocol adds an arbitrary
offset value and/or the result of its iteration, in an effort to avoid
revealing its own initial value or the initial values of other nodes. What is
important is for each node to ensure that the total offset it adds cancels out
in the end (the accumulated sum of offsets is zero). The work in [13]
establishes that under certain conditions on the communication topology, the
protocol allows the nodes to calculate the average of their initial values in
a privacy-preserving manner, despite the presence of curious agents. For
directed topologies, inspired by ratio consensus strategies, we devised in
[22] a distributed mechanism, where each node updates its information state by
combining the available information received from its in-neighbors using
constant positive weights and by adding an offset to one of the two states
communicated during the execution of the algorithm. As shown in [22], if after
a finite number of iteration the nodes ensure that the total offset is zero,
this privacy-preserving version of ratio consensus converges to the exact
average of the nodes’ initial values, even in the presence of bounded time-
varying delays.

The privacy-preserving gossip-based protocol we develop and analyze in this
paper enables all of the nodes to calculate in finite time an $\epsilon$-close
approximation of the average of their initial values, without loss of privacy
despite the presence of multiple (possibly colluding) curious nodes. Curious
nodes are assumed to have full knowledge of the proposed protocol and are
allowed to collaborate arbitrarily among themselves, but do _not_ interfere in
the computation of the average value of the network in any other way. Our
approach does not depend on any cryptographic algorithm, but operates by
allowing the nodes to introduce pseudo-random offsets (unknown to the curious
nodes). Specifically, the proposed protocol enhances the standard gossip-based
protocol [9] in two main directions. The first enhancement is that, at each
time-step, the randomly chosen pair of nodes not only exchanges information
and calculates the local (pairwise) average value, but if any of the two is
nodes following the protocol, then that node also adds an arbitrary offset
value to its updated (average) result, in an effort to avoid revealing its own
initial value or the initial values of other nodes. What is important is for
each node to ensure that the total offset added (accumulated sum of offsets)
eventually becomes zero. We establish simple conditions to ensure that the
proposed protocol allows the nodes to calculate the average of their initial
values in a privacy-preserving manner, despite the presence of curious agents.
In particular, the protocol guarantees that when there are at least two nodes
following the protocol then privacy is guaranteed for both of them in the
sense that their individual initial values are not revealed to the curious
nodes. Note that it might still be possible for the curious nodes to determine
exactly the sum of the initial values of nodes that follow the privacy-
preserving protocol (but not their individual values). The second enhancement
is a deterministic criterion that enables the nodes to determine, in a
distributed manner, when to seize initiating gossip exchanges because
approximate average consensus has been reached, i.e., it is shown that
eventually all nodes will seize activity, and that point they will all have
values that are within a small distance from the true average of their initial
values.

The remainder of the paper is organized as follows. In Section II, we provide
some relevant background on graph theory and describe the linear and gossip-
based iterative strategies for reaching average consensus in a given
distributed system. In Section III we formulate the problem of interest. We
describe our proposed privacy-preserving gossip protocol, and the main results
of the paper in Section IV. In Section V we present an example, and finish the
paper with conclusions and directions for future work in Section VI.

##  II Background

###  II-A Distributed System Model

In a distributed system we can model the network topology as a directed graph
(digraph) G={X,E} where X={$1,2,...,N$} is the set of components in the system
and $E\subseteq X\times X-\\{(i,i)\mid{i\in X}$} is the set of edges (self-
edges excluded). In particular, edge $(i,j)\in E$ if node $j$ can send
information to node $i$. In this work, we focus on components that interact
via bidirectional links in a way that forms a connected undirected (or
symmetric) graph, i.e., a graph for which $(i,j)\in E$ if and only if
$(j,i)\in E$. The set of nodes that can exchange information with node $i$ are
said to be the neighbors of node $i$ and are represented by the set
$\mathcal{N}_{i}=\\{(j,i)\in E\\}$. The number of neighbors of node i is
called the degree of node $i$ and is denoted as
$\mathcal{D}_{i}$=$|\mathcal{N}_{i}|$.

Our model deals with networks where information is transmitted/received
between a randomly selected pair of nodes at each time step. More
specifically, each node randomly wakes up and initiates an information
exchange with a randomly selected neighbor. We assume that nodes are aware of
the unique identity (id) of each node they communicate with and that, during
the exchange process, all the information to/from a particular node is
transmitted/received successfully [3]. Moreover, the nodes must have
sufficient memory and computational capability in order to store and perform
certain computations while the iteration is executing. During the
transmission/reception process, the nodes in the network receive a value from
the selected node, and transmit their value back to the selected node.

###  II-B Average Consensus via Linear Iterative Strategy

In average consensus problems the objective is the calculation of the average
of the initial values of the nodes in the network. Assume that each node $i$
has some initial value $x_{i}$ and, at each time-step $k$, it updates its
value as a weighted sum of its own value and the values of its in-neighbors
(e.g., following the method in [7]). Specifically, at each time-step $k$, each
node updates its value as

$\displaystyle
x_{i}[k+1]=w_{ii}x_{i}[k]+\sum_{j\in\mathcal{N}_{i}}w_{ij}x_{j}[k],$ (1)

where $w_{ij}$ are a set of (fixed) weights and $x_{i}[0]=x_{i},\forall i\in
X$. The values for all the nodes at time-step $k$ can be aggregated into the
value vector $x[k]=\big{[}x_{1}[k]$ $x_{2}[k]$ $...$ $x_{N}[k]\big{]}^{T}$ and
the update strategy for the entire network can be written compactly as

$x[k+1]=\underbrace{\left[\begin{array}[]{cccc}w_{11}&w_{12}&\cdots&w_{1N}\\\
w_{21}&w_{22}&\cdots&w_{2N}\\\ \vdots&\vdots&\ddots&\vdots\\\
w_{N1}&w_{N2}&\cdots&w_{NN}\end{array}\right]}_{W}x[k],$

for $k\in\mathbb{N}$, where $w_{ij}=0$ if
$x_{j}\notin\textit{$\mathcal{N}_{i}$}\cup{\\{i\\}}$.

The system is said to reach asymptotic consensus if
$lim_{k\rightarrow\infty}x_{i}[k]=\textit{f}(x_{1}[0],x_{2}[0],...,x_{N}[0])$
for each i, where $\textit{f}(x_{1}[0],x_{2}[0],...,x_{N}[0])\in\mathbb{R}$.
When $\textit{f}(x_{1}[0],x_{2}[0],...,x_{N}[0])=c^{T}x[0]$ for some column
vector $c$ (where $c^{T}$ is the transpose of $c$), the following result by
Xiao and Boyd from [23] characterizes the conditions under which iteration (1)
achieves asymptotic consensus.

###### Theorem 1

[23] The iteration given by (1) reaches asymptotic consensus on the linear
functional $c^{T}x[0]$ (under the technical condition that $c$ is normalized
so that $c^{T}\bf{1}=1$, where $\bf{1}$ is the all ones column vector of size
$N$) if and only if the weight matrix $W$ satisfies the conditions below:

  1. 1.

$W$ has a simple eigenvalue at 1, with left eigenvector $c^{T}$ and right
eigenvector $\bf{1}$=[1 1…1$]^{T}$;

  2. 2.

All other eigenvalues of $W$ have magnitude strictly less than 1.

In particular, if $c^{T}$=$\frac{1}{N}$[1 1…1], then average consensus is
reached. Also note that if $w_{ij}$ are restricted to be nonnegative, then the
above conditions are equivalent to $W$ being a doubly stochastic matrix.

###  II-C Average Consensus using a Randomized Gossip Strategy

In a randomized gossip protocol, a pair of nodes is randomly selected at each
iteration and perform local (pairwise) averaging of the two values.
Specifically, at iteration $k$, nodes $i$ and $j$ (if selected) update their
values as

$\displaystyle x_{i}[k+1]=x_{j}[k+1]=\frac{x_{i}[k]+x_{j}[k]}{2},$ (2)

with $x_{i}[0]=x_{i}$ for all $i$. Typically, it is assumed that each node
“wakes up” in an asynchronous manner and randomly selects a neighbor to
perform local averaging with. Note that the remaining nodes do not update
their values, i.e., $x_{l}[k+1]=x_{l}[k]$ for all nodes $l$, $l\neq j$ and
$l\neq i$.

It can be shown that, if nodes are following the above protocol, they are
guaranteed to asymptotically reach the average of their initial values
$\bar{x}=\frac{1}{N}\sum_{l=1}^{N}x_{l}[0]$, as long as the network is
connected (so that all the available information can flow to every node in the
network) and all pairs of nodes in the network are active frequently enough.

###  II-D Previous Work on Privacy-Preserving Average Consensus

Privacy-preserving average consensus in the presence of curious agents in the
network has already received significant attention in the literature (e.g.,
[10], [13], [9], [14] [12]). Below, we briefly review the most relevant
existing results for our work here.

As mentioned earlier, the authors of [10] proposed a transformation method
using random offset values in a cooperative wireless network. Specifically,
each node $i$ that wishes to protect its privacy adds a random offset value
$u_{i}$ to its initial value $x_{i}$. This ensures that its value will not be
revealed to curious nodes that might be observing the exchange of values in
the network. The idea is based upon the observation that, when a large number
of nodes employ the protocol, their offsets will have a zero net effect on the
average with high probability, allowing the nodes to converge to the true
average of the initial values of the nodes in the network. Specifically, each
node $i$ sets $x_{i}[0]$=$x_{i}$′= $x_{i}$+$u_{i}$ where $u_{i}$,
$i=1,2,...,N,$ are i.i.d. random variables with zero mean. Then, following the
protocol for asymptotic average consensus, the nodes converge to,

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}x_{i}^{\prime}=\underbrace{\frac{1}{N}\sum_{i=1}^{N}x_{i}}_{\overline{x}}+\underbrace{\frac{1}{N}\sum_{i=1}^{N}u_{i}}_{U},$
(3)

where $\overline{x}$ is the desirable average of the original initial values
and $U$ is a random variable that captures the net effect of the offsets.
Clearly, $E[U]=0$ (since the $u_{i}$ are zero mean) and
$var[U]=\frac{1}{N}var(U_{i})$ (since the $u_{i}$ are i.i.d.).

For ${N\to\infty}$ we have $var(U)~{}\rightarrow 0$ which means that

$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}x^{\prime}_{i}=\lim_{N\to\infty}\bigg{[}\frac{1}{N}\sum_{i=1}^{N}x_{i}+\frac{1}{N}\sum_{i=1}^{N}u_{i}\bigg{]}=\overline{x}+0=\overline{x},$
(4)

where the last inequality is taken in the mean square sense. For a large
number of nodes $(N\rightarrow\infty)$, this method can give results very
close to the true average of the network; however, as the number of nodes
decreases, the accuracy of this method also decreases, due to the fact that
the offset values added to the protected nodes will add a random offset (with
mean zero and some finite variance) to the true average value of the network.

In our own previous work [13], the objective is also to calculate the average
of the initial values of the nodes in the network, while at the same time
preserving the privacy of the nodes following the protocol. The scheme that is
used in [13] assumes that the underlying network forms a connected undirected
graph and makes use of linear iterations as in (1) were the weights $w_{ij}$
form a doubly stochastic matrix $W=[w_{ij}]$ (thus, the nodes asymptotically
reach consensus to the average of their initial values). The main difference
is that node $i$ following the protocol sets its initial value
$x^{\prime}_{i}[0]=x_{i}+u_{i}$ (where $u_{i}$ is some random offset), and
subsequently updates its value as

$\displaystyle
x_{i}^{\prime}[k+1]=w_{ii}[k]x_{i}^{\prime}[k]+\sum_{j\in\mathcal{N}_{i}}w_{ij}x_{j}^{\prime}[k]+u_{i}[k],\
k=0,1,...,$ (5)

where $u_{i}[k]$ is a pseudo-random value chosen by node $i$ at time-step $k$.
The constraint is that $u_{i}[k]=0$ for $k>L_{i}$ (for some $L_{i}$ known only
to node $i$) and

$\displaystyle u_{i}[L_{i}]=-\sum_{t=0}^{L_{i}-1}u_{i}[t]-u_{i}.$ (6)

At time-step $L_{i}$, node $i$ effectively cancels-out the pseudo-random
values it has added during the information exchange in the network up to that
point.

The main contribution of [13] is the establishment of topological conditions
that ensure privacy for the nodes following the proposed protocol despite the
presence of curious nodes in the network. Considering a fixed network with $N$
nodes described by a digraph $G=\\{X,E\\}$ and the iteration in (1) with
weights that form a primitive doubly stochastic weight matrix $W$, it assumes
that a set of nodes $P$ follow the predefined privacy-preserving strategy in
(2) with random offsets chosen as in (6). Curious node $i$ will not be able to
identify the initial value $x_{j}$ for $j\in{P}$, as long as $j$ has at least
one other node $l$ connected to it for which all paths from $l$ to the curious
node $i$ are through a node $j^{\prime}$ following the protocol (i.e.,
$j^{\prime}\in P$). Specifically, if the above mentioned condition is
satisfied, the network will reach average consensus and the privacy of the
initial values of the nodes following the protocol will be preserved during
the linear iteration process.

##  III Problem Setup

Consider a network with $N$ nodes each with some initial value (node $i$ has
initial value $x_{i}$). Let us assume that all nodes follow a gossip-based
strategy for asymptotically reaching agreement to the average of their initial
values. However, we consider that (i) a set of curious nodes $X_{c}\subset X$
try to identify the initial values of all or some of the nodes in the network,
(ii) a set of private nodes $X_{p}\subset X$ would like to preserve their
privacy by not revealing to other nodes their initial values, and (iii) a set
of neutral nodes $X_{n}\subset X$ are neither curious nor try to preserve the
privacy of their initial values (we allow some nodes not to be curious and not
follow the privacy-preserving protocol in order to investigate the worst-case
scenario that this protocol can handle). In this paper we describe a
distributed protocol that allows each node $i$ to calculate the average
${\overline{x}}$, while every node $i\in X_{p}$ is able to preserve the
privacy of its initial value from possibly colluding curious nodes in $X_{c}$.
Furthermore, once values close to $\overline{x}$ obtained for all nodes, each
node is able to determine when to terminate its operation in a distributed
manner using a deterministic criterion. More specifically when the nodes
terminate their operation, we are guaranteed that approximate average
consensus has been reached, i.e., all nodes have obtained values that are
within a small distance from $\overline{x}$. Note that in this paper, we
assume that the curious nodes have full knowledge of the proposed protocol and
are allowed to collaborate arbitrarily among themselves (exchanging
information as necessary), but do not interfere in the computation of the
average value in any other way.

##  IV Proposed Strategy and Main Results

The main contribution of this paper includes a deterministic privacy-
preserving gossip protocol that also guarantees distributed stopping. The
protocol is deterministic in the sense that, at its completion, it is
guaranteed that all values are close to $\overline{x}$ and the nodes that
follow the privacy-preserving protocol have protected their initial values
from exposure, as long as a simple sufficient condition is satisfied.

###  IV-A Deterministic Privacy-Preserving Protocol

The first objective of the protocol is to calculate the average of the initial
values of the nodes in the network while preserving the privacy of the nodes
that follow the proposed protocol. The scheme that we study makes use of the
iteration in (2) where, at each iteration, a randomly selected pair of nodes
perform local averaging of their values and asymptotically reach approximate
consensus to the average of their initial values $\overline{x}$, at least as
long as each pair of nodes is selected enough times. What we propose in this
work is a deterministic privacy-preserving protocol in which node $i$
following the protocol sets its initial value to

$\displaystyle x^{\prime}_{i}[0]=x_{i}+u_{i}$ (7)

where $u_{i}$ is some random offset. A node $l$ that is not following the
protocol simply sets $x^{\prime}_{l}[0]=x_{l}$. Subsequently, when nodes $i$
and $j$ are selected to perform local averaging (say at iteration $k$), node
$i$ updates its value as

$\displaystyle
x_{i}^{\prime}[k+1]=\big{(}x_{i}^{\prime}[k]+x_{j}^{\prime}[k]\big{)}/2+u_{i}[k],\
k=0,1,...,$ (8)

where $u_{i}[k]$ is a pseudo-random value chosen by node $i$ at iteration $k$.
Node $j$ does something similar (adding a pseudo-random offset $u_{j}[k]$) if
it is following the protocol; if node $j$ is not following the protocol, then
node $j$ sets

$\displaystyle
x_{j}^{\prime}[k+1]=\big{(}x_{i}^{\prime}[k]+x_{j}^{\prime}[k]\big{)}/2,\
k=0,1,....$ (9)

The constraints are (i) $u_{i}[k]=0$ for $k>L_{i}$ where $L_{i}$ is the time
step where node $i$ is first selected after it identifies that it has
exchanged information with each node in the network at least once, and (ii)

$\displaystyle u_{i}[L_{i}]=-\sum_{t=0}^{L_{i}-1}u_{i}[t]-u_{i}.$ (10)

Note that at time-step $L_{i}$, node $i$ effectively cancels-out the pseudo-
random values it has added during the information exchange in the network up
to that point; thus, from that point onward the remaining accumulated offset
(due to node $i^{\prime}s$ offset values) in the system is zero. Note that for
notational convenience we take $u_{i}[k]=0$ for all time steps $k$ at which
node $i$ is not selected to update its value.

Protocol Description: We now describe the distributed protocol from the
perspective of node $i$. When selected, nodes following the protocol run the
linear iteration in (8) or (9) (depending on whether they follow the protocol
or not). Each time a link is activated, each node among the selected pair of
nodes is aware of the identity of the other node involved; thus, nodes are
able to identify when they have communicated at least once with all the nodes
in the network. We refer to $L_{i}$ as the iteration ($k=L_{i}$) at which node
$i$ is first selected after it has communicated at least once with all other
nodes. A node not following the protocol sets $u_{i}=0$ and $u_{i}[k]=0$ for
$k=0,1,2,...,$ which is the standard gossip protocol for reaching average
consensus. Node $i$ that follows the privacy-preserving protocol executes (8),
with $x^{\prime}_{i}[0]=x_{i}+u_{i}$ and

  * i)

Chooses a pseudo random offset $u_{i}[k],\ k=0,1,...,L_{i}-1$ for integer
$L_{i}$, ($u_{i}[t]=0$ for iteration steps $t$ for which node $i$ is not
selected);

  * ii)

Sets

$\displaystyle u_{i}[L_{i}]=-\sum_{t=0}^{L{{}_{i}}-1}u_{i}[t]-u_{i};$ (11)

  * iii)

Sets $u_{i}[k]=0$ for $k>L_{i}$.

Note that $L_{i}$ is an integer number of steps that becomes known to each
node i following the protocol. This number is a random variable since the
pairs are selected randomly at each iteration.

$Lemma~{}1$: Following the iteration in (1) and in combination with the
constraint in (11) the network will reach asymptotic average consensus.

###### Proof:

It is not hard to see that, if we let $L_{max}=max_{i}\\{L_{i}\\}$, then

$\displaystyle\sum_{i=1}^{N}x^{\prime}_{i}[L_{max}+1]=\sum_{i=1}^{N}x_{i};$

then, using the fact that from this point onward nodes follow a standard
gossip protocol, we conclude that

$\displaystyle\lim_{k\to\infty}x^{\prime}_{i}[k]=\frac{1}{N}{\sum_{i=1}^{N}x^{\prime}_{i}[L_{max}+1]}=\frac{1}{N}{\sum_{i=1}^{N}x_{i}}=\overline{x},$

for all $i\in X$ (i.e., average consensus has been reached).∎

1: Consider a fully connected graph $G=\\{X,E\\}$ with $N=|X|$ nodes. Each
node $i$ has initial value as $x_{i}$.
2: If node $i$ is following the protocol, it sets its initial value to
$x^{\prime}_{i}[0]=x_{i}+u_{i}$ (where $u_{i}$ is some random offset);
otherwise, $x_{i}^{\prime}[0]=x_{i}$.
3: Let $L_{i}$ be the time step at which node $i$ is next activated after it
communicates at least once with all the other nodes in the network). During
the first $L_{i}-1$ steps, if node $i$ is selected, it adds a random offset
value to its transmitted value. At time-step $L_{i}$, node $i$ cancels out the
injected noise in the system by setting
$u_{i}[L_{i}]=-\sum_{t=0}^{L{{}_{i}-1}}u_{i}[t]-u_{i}$, and for $k>L_{i},$ it
sets $u_{i}[k]=0$.
4: Specifically, for $k=0,1,2,...,$ each node $i\in X$, if selected with node
$j$, does the following ($u_{i}[k]=0$ for all $k$ if $i$ is not following the
protocol)
5: if $k<L_{i}$ then:
$x_{i}^{\prime}[k+1]=\frac{(x_{i}^{\prime}[k]+x_{j}^{\prime}[k])}{2}+u_{i}[k]$
6: if $k=L_{i}$ then:
$x_{i}^{\prime}[k+1]=\frac{(x_{i}^{\prime}[k]+x_{j}^{\prime}[k])}{2}-\sum_{t=0}^{L{{}_{i}-1}}u_{i}[t]-u_{i}$
7: if $k>L_{i}$ then:
$x_{i}^{\prime}[k+1]=\frac{x_{i}^{\prime}[k]+x_{j}^{\prime}[k]}{2}$
8: End

Algorithm 1 Privacy-Preserving Gossip Protocol

###  IV-B Deterministic Distributed Stopping Protocol

In the proposed deterministic distributed stopping algorithm, each node makes
a decision on how to update and/or transmit its value, based on the difference
between its calculated value and the value it receives from its selected
neighbor at each iteration that the node is activated. In the algorithm, the
iterative process ends when all nodes cease transmitting their values, in
which case they can be shown to have reached approximate average consensus,
i.e., the absolute difference between the final value of each node and the
exact average $\overline{x}$ of the initial values is smaller than an error
bound $\epsilon$ (small real value).
In terms of the definition below we are interested in reaching
$\varepsilon$-approximate average consensus and also in identifying (in a
distributed manner) when such approximate average consensus has been reached.

$Definition~{}2~{}(\varepsilon-Approximate~{}Average~{}Consensus):$ At the end
of the iterative process, the nodes have reached $\varepsilon$-approximate
average consensus if the value $x_{i}[f]$ of each node $i\in X$ satisfies
$|x_{i}[f]-\bar{x}|\leq\varepsilon\;,\;\forall i\in X$, where $\bar{x}$ is the
average of the initial values.

Protocol Description: We now describe the distributed stopping protocol from
the perspective of node $i$. Nodes following the protocol run the linear
iteration in (8) and (9), depending on whether or not they run the privacy-
preserving protocol, in order to reach average consensus. Nodes running the
privacy-preserving protocol wait until they communicate at least once with
each other node and then start executing the distributed stopping protocol
below. The protocol has two phases. During phase 1, each time node $i$ is
activated (say with node $j$) checks whether the absolute difference
$|x_{i}[k]-x_{j}[k]|<\varepsilon$; if the condition is satisfied then the
flags $f_{ji}$ and $f_{ij}$ that represent the closeness of the values between
the two nodes become $1$; if the condition is not satisfied then the nodes
proceed to perform local averaging of their values and set all their flags to
zero. Once node $i$ has all its flags $\\{f_{ji}\mid j\in X\\}-\\{f_{ii}\\}$
equal to $1$, it enters the second phase. In this phase, node $i$ does not
initiate any exchanges but, if it is selected by another node $j$, it checks
whether the absolute difference $|x_{i}[k]-x_{j}[k]|>\varepsilon$ holds; in
such case, all the flags of node $i$ become $0$ and node $i$ goes back to
phase 1.

Algorithm 2 Distributed Stopping for Average Consensus

1:Input: Each node $i$ initializes $x_{i}[0]=x_{i}$ and if node $i$ follows
the privacy protocol, then it starts to implement the distributed stopping
protocol described below, after time step $L_{i}$.

2:Initialize all flags $\\{f_{ji}\mid j\in X-\\{i\\}\\}$ to zero.

3:For each $k>L_{i}$ when a pair of nodes $i$ and $j$ is randomly selected,
nodes $i$ and $j$ do the following.

4: If $\mid x_{i}[k]-x_{j}[k]\mid<\varepsilon$

5: $Flags_{ij}=Flags_{ji}=1$

6: else

7: $x_{i}[k+1]=x_{j}[k+1]=\frac{x_{i}[k]+x_{j}[k]}{2}$

8: $Flags_{*i}=Flags_{*j}=0$ (set all flags to zero)

9: End

10: If all $f_{ji}=1$ node $i$ does not initiate any gossip interactions
(unless contacted by another node $j$).

11: End

12:End

###### Theorem 2

Consider a network described by a fully connected undirected graph
$G=\\{X,E\\}$, where nodes run the enhanced gossip algorithm in order to reach
average consensus. Following the deterministic distributed stopping protocol,
nodes reach $\varepsilon$-approximate average consensus after a finite number
of iterations.

###### Proof:

Due to space limitations, we only provide a sketch of the proof. We can first
establish, by contradiction, that the nodes will stop after a finite number of
iterations $f$. Suppose that the iteration runs forever; this means that at
each iteration $k$ at least two neighboring nodes are active (transmitting
nodes). Let
$A=\\{(i_{1},j_{1}),(j_{1},i_{1}),(i_{2},j_{2}),(j_{2},i_{2}),\ldots,(i_{\kappa},j_{\kappa}),(j_{\kappa},i_{\kappa})\\}\subseteq
E$ be the non-empty set of edges that are active infinitely often. Also, let
$\ell$ be the latest iteration at which an edge in the set $E-A$ is active
(note that $\ell$ is finite). We have that the active edge at iteration $k>l$
denoted by $e[k]$ satisfies $e[k]\in A$ (i.e., after time step $\ell$ each
active edge belongs in the set $A$). Since each $e[k]$ results in an
undirected graph $G[k]=\\{X,\\{e[k\\}\\}$, the graph $\\{X,A\\}$ is also an
undirected graph and can be partitioned into connected components. Then, from
[15], it follows that the nodes in each connected component asymptotically
reach the same value. This implies that the edges in the connected component
would cease to be active, which is a contradiction. Since the nodes stop at
some iteration $f$, we have $\varepsilon$-approximate local consensus [24].
Using Proposition 2 from [24], we also have $(D\times\varepsilon)$-approximate
global consensus, and since we have a fully connected graph $(D=1)$, the
absolute difference between any pair of node values after they stop
transmitting is at most $\varepsilon$. Moreover, we have

$x_{\min}[f]\leq\frac{1}{N}\sum_{l=1}^{N}x_{l}[f]\leq x_{\max}[f]$

(where $x_{\min}[f]=\min_{i}\left\\{x_{i}[f]\right\\}$ and
$x_{\max}[f]=\max_{i}\left\\{x_{i}[f]\right\\}$), which implies that

$x_{\min}[f]\leq\bar{x}=\frac{1}{N}\sum_{l=1}^{N}x_{l}=\frac{1}{N}\sum_{l=1}^{N}x_{l}[f]\leq
x_{\max}[f]\;.$

Since $\mid x_{max}[f]-x_{min}[f]\mid<\varepsilon$ this completes the proof. ∎

###  IV-C Conditions for Privacy-Preservation

###### Theorem 3

Let us consider a fully connected graph $G=\\{X,E\\}$ of $N$ nodes. The set of
private nodes $X_{p}$ follow the predefined privacy-preserving gossip strategy
in (2) with random offsets chosen as in (8) and (11). The set of curious nodes
$X_{c}$ will not be able to identify the initial value $x_{j}$ of node $j\in
X_{p}$, as long as

  1. 1.

there is at least one node $l\in X_{p}$ in the network that also follows the
proposed privacy-preserving protocol;

  2. 2.

there is at least one node $l\in X_{n}$ in the network which first exchanges
information with node $j$ (at the first iteration it is involved in a selected
pair).

Note that, if the conditions in Theorem $2$ are also satisfied, the network
will reach approximate average consensus in finite time, while the privacy of
the initial values of the nodes following the protocol will be preserved
during the process.

###### Proof:

The proof of Theorem 3 analyzes the following simple scenarios.

  * 1)

If all nodes $i\in X-\\{j\\}$ in the network are curious and they communicate
with each other, it is not possible for $j$ to maintain the privacy of its
initial value $x_{j}$. Specifically, at the initialization of our algorithm,
the curious nodes will know the value $x_{j}+u_{j}$ (recall that $u_{j}$ is
the value that $j$ injected to its initial value in (7)). During the iteration
of our algorithm, curious nodes in $X_{c}$ will get to learn the values
$u_{j}[k]$ for $k=0,1,...,L_{j}$, that $j$ injects to $x_{j}^{\prime}[k]$
(since curious nodes collude with each other and exchange information, they
can compare each pair of consecutive values node $j$ announces to determine
$u_{j}[0]$, then $u_{j}[1]$, and so forth up to $u_{j}[L_{j}]$). This means
that the curious nodes $i\in X_{c}$ will be able to compute the total amount
of injected offsets $u_{j}+\sum_{t=0}^{L_{j}-1}u_{j}[t]=-u_{j}[L_{j}]$ and,
from that, $u_{j}$ and therefore the initial value $x_{j}$ of node $j$. As a
result, the privacy of node $j$ is not preserved and at least one node (other
than $j$) that is not curious is needed in the network. It is also worth
highlighting that no privacy-preserving protocol can protect node $j$ in this
case (since $\overline{x}$ becomes known and the curious nodes know their own
value).

  * 2)

Let us consider the case for which there exists at least one node $l\in X_{n}$
in the network (i.e., node $l$ is neither curious nor following the privacy-
preserving protocol). Also, $i\in X_{c},\ \forall\ i\in X\setminus\\{l,j\\}$
(i.e., all other nodes in the network are curious). At the initialization of
our algorithm, node $j$ will inject an initial offset to its value
$x_{j}^{\prime}[0]=x_{j}+u_{j}$ as in (7) while $l$ will maintain its value
$x_{l}$, i.e., $x_{l}^{\prime}[0]=x_{l}$. Let us assume that at the first time
step $k=0$, $j$ and $l$ are the selected pair. Then, they will update their
values as

$\displaystyle x_{j}^{\prime}[1]=$
$\displaystyle\big{(}x_{j}^{\prime}[0]+x_{l}^{\prime}[0]\big{)}/2+u_{j}[0],\ $
and $\displaystyle x_{l}^{\prime}[1]=$
$\displaystyle\big{(}x_{j}^{\prime}[0]+x_{l}^{\prime}[0]\big{)}/2.$ (12)

When nodes $j$ and $l$ are paired with curious nodes at time steps
$k=1,2,...$, their values $x_{j}^{\prime}[1]$ and $x_{l}^{\prime}[1]$ will
become known (along with $(x_{j}^{\prime}[0]+x_{l}^{\prime}[0])/2$ and
$u_{j}[0]$). Later on, curious nodes could determine $u_{j}[1]$, $u_{j}[2]$
and eventually $u_{j}[L_{j}]=-u_{j}-\sum_{t=1}^{L_{j}-1}u_{j}[t]$. Thus, they
can determine $u_{j}$ and compute the exact sum of the initial values of nodes
$j$ and $l$ as $x_{j}+x_{i}=2x_{l}^{\prime}[1]-u_{j}$ (the value of
$x_{l}^{\prime}[1]$ becomes known at iteration $k=1$, whereas $u_{j}$ can be
calculated after iteration $L_{j}$, once $u_{j}[L_{j}]$ is available).

  * 3)

Let us consider the case for which there exists at least one node $l\in X_{p}$
in the network that also follows the proposed privacy-preserving protocol.
Also, $i\in X_{c},\ \forall\ \in X\setminus\\{j,l\\}$ (i.e., all other nodes
in the network are curious). At the initialization of our algorithm, $j$ and
$l$ will inject initial offsets to their values
$x_{j}^{\prime}[0]=x_{j}+u_{j}$ and $x_{l}^{\prime}[0]=x_{l}+u_{l}$,
respectively, as in (7). During the iteration of our algorithm, each time $j$
or $l$ are paired with another curious node $i\in X_{c}$ at time steps
$k=1,2,...$, they will inject offsets to their values $u_{j}[k]$, and
$u_{l}[k]$ as

$\displaystyle x_{j}^{\prime}[k+1]=$
$\displaystyle\big{(}x_{j}^{\prime}[k]+x_{i}^{\prime}[k]\big{)}/2+u_{j}[k],\ $
or $\displaystyle x_{l}^{\prime}[k+1]=$
$\displaystyle\big{(}x_{l}^{\prime}[k]+x_{i}^{\prime}[k]\big{)}/2+u_{l}[k].$
(13)

Let us assume that at time step $k^{\prime}$, nodes $j$ and $l$ are selected
as a pair. Both nodes will inject offsets $u_{j}[k^{\prime}]$ and
$u_{l}[k^{\prime}]$ to their values as

$\displaystyle x_{j}^{\prime}[k^{\prime}+1]=$
$\displaystyle\big{(}x_{j}^{\prime}[k^{\prime}]+x_{l}^{\prime}[k^{\prime}]\big{)}/2+u_{j}[k^{\prime}],\
$ and $\displaystyle x_{l}^{\prime}[k^{\prime}+1]=$
$\displaystyle\big{(}x_{l}^{\prime}[k^{\prime}]+x_{j}^{\prime}[k^{\prime}]\big{)}/2+u_{l}[k^{\prime}].$
(14)

The sum of values of $u_{j}[k^{\prime}-1]$ and $u_{l}[k^{\prime}-1]$, which
were injected after the previous exchange of each node and are known only to
nodes $j$ and $l$, additively affect both $x_{j}^{\prime}[k^{\prime}+1]$ and
$x_{i}^{\prime}[k^{\prime}+1]$. As a result, the privacy of both nodes $j$ and
$l$ is preserved (curious nodes can only resolve the sum of
$u_{j}[k^{\prime}-1]+u_{l}[k^{\prime}-1]$, but not the individual values). As
we mentioned in the previous scenario, since only nodes $j$ and $l$ are
following the proposed protocol, curious nodes $i\in X_{p}$ will be able to
compute the exact sum of the initial values $x_{j}+x_{l}$ of nodes $j$ and
$l$.

From the above scenarios, we have that node $j$ is able to preserve the
privacy of its initial value if (i) there exists at least one node $l\in
X_{n}$ in the network and node $l$ first exchanges information with node $j$
among them, or (ii) there exists at least one other node $l\in X_{p}$, $l\neq
j$, that also follows the proposed privacy-preserving protocol. In both cases,
the privacy of the initial value of node $j\in X_{p}$ is preserved.
This completes the proof. ∎

##  V Illustrative Example

In this section, we present an example to illustrate the behavior and
potential advantages of our proposed algorithm. We analyze the case of a fully
connected graph of $20$ nodes with the initial value $x_{i}$ for each of the
nodes chosen independently and uniformly in $[0,1]$. For this case, each node
$i\in X$ executes Algorithm 1 and Algorithm 2, i.e., it utilizes (i) the
privacy-preserving protocol described in Section IV-A, and (ii) the
distributed stopping protocol described in Section IV-B. Specifically, during
the operation of Algorithm 1, in order to ensure privacy-preservation, each
node $i$ infuses pseudo random offsets $u_{i}[k]$ chosen uniformly in
$[-1,1]$, for $\ k=0,1,...,L_{i}-2$ where $L_{i}$ is the number of steps
elapsed at the first activation of node $i$ after it has communicated at least
once with all the other nodes in the network (note that $u_{j}[k]$ is zero
when node $j$ is not selected). For distributed stopping guarantees, each node
$i$ initially sets its flag to $0$ and takes the value of
$\varepsilon=0.0001$.

In Fig. 1, the average of the initial states of the nodes is equal to
$\bar{x}=0.3699$. We show the value of each node $i$ plotted against the
required time steps for convergence. We observe that each node calculates the
approximate average of the initial states after $2140$ time steps. When the
nodes stop, the error bound between the absolute difference of the value of
any node in the system from the exact average value ${\overline{x}}$ is
smaller than $\varepsilon$.

Figure 1: Convergence of gossip-based privacy-preserving protocol in finite
time.

##  VI CONCLUSIONS AND FUTURE WORK

In this paper, we have considered the problem of privacy-preserving average
consensus in finite time in the presence of curious nodes. We presented a
distributed privacy-preserving average consensus algorithm that enables all of
the components of a distributed system, each with some initial value, to reach
(approximate) average consensus on their initial values after executing a
finite number of iterations. We proposed an enhancement of the standard gossip
protocol that can be followed by each node that does not want to reveal its
initial value to curious nodes with which it might interact. We characterized
the conditions on the information exchange that guarantee privacy-preservation
for a specific node. Furthermore, we also provided a criterion that allows the
nodes to determine, in a distributed manner, when to terminate their operation
because approximate average consensus has been reached, Finally, we presented
simulation results of our proposed protocol.

In the future we plan to extend the results of this paper to the case where
the proposed protocol operates over arbitrary directed graphs (not necessarily
fully connected). Furthermore, achieving distributed stopping for finite time
average consensus with privacy guarantees is to be investigated in our future
work.

## References

  * [1] A. D. Domínguez-García and C. N. Hadjicostis, “Distributed algorithms for control of demand response and distributed energy resources,” in _Proceedings of $50^{th}$ IEEE Conference on Decision and Control and European Control Conference_, 2011, pp. 27–32.
  * [2] W. Ren, R. W. Beard, and E. M. Atkins, “A survey of consensus problems in multi-agent coordination,” in _Proceedings of American Control Conference_ , 2005, pp. 1859–1864.
  * [3] N. Lynch, _Distributed Algorithms_. San Mateo: CA: Morgan Kaufmann Publishers, 1996.
  * [4] R. Koetter and M. Médard, “An algebraic approach to network coding,” _IEEE/ACM Transactions on Networking_ , vol. 11, no. 5, pp. 782–795, October 2003.
  * [5] J. Hromkovic, R. Klasing, A. Pelc, P. Ruzicka, and W. Unger, _Dissemination of Information in Communication Networks_ , Springer-Verlag, 2005.
  * [6] J. Cortés, “Distributed algorithms for reaching consensus on general functions,” _Automatica_ , vol. 44, no. 3, pp. 726–737, March 2008.
  * [7] S. Sundaram and C. N. Hadjicostis, “Distributed function calculation and consensus using linear iterative strategies,” _IEEE Journal on Selected Areas in Communications_ , vol. 26, no. 4, pp. 650–660, May 2008.
  * [8] C. N. Hadjicostis, A. D. Domínguez-García, and T. Charalambous, “Distributed averaging and balancing in network systems, with applications to coordination and control,” _Foundations and Trends® in Systems and Control_ , vol. 5, no. 3–4, 2018.
  * [9] A. G. Dimakis, S. Kar, J. M. F. Moura, M. G. Rabbat, and A. Scaglione, “Gossip algorithms for distributed signal processing,” _Proceedings of the IEEE_ , vol. 98, no. 11, pp. 1847–1864, November 2010.
  * [10] M. Kefayati, M. S. Talebi, B. H. Khalaj, and H. R. Rabiee, “Secure consensus averaging in sensor networks using random offsets,” in _IEEE International Conference on Telecommunications and Malaysia International Conference on Communications_ , May 2007, pp. 556–560.
  * [11] E. Nozari, P. Tallapragada, and J. Cortes, “Differentially private average consensus: Obstructions, trade-offs, and optimal algorithm design,” _Automatica_ , vol. 81, pp. 221–231, 2017.
  * [12] J. Cortés, G. E. Dullerud, S. Han, J. L. Ny, S. Mitra, and G. J. Pappas, “Differential privacy in control and network systems,” in _$55^{th}$ IEEE Conference on Decision and Control_, 2016, pp. 4252–4272.
  * [13] N. Manitara and C. N. Hadjicostis, “Privacy-preserving asymptotic average consensus,” in _European Control Conference_ , Zurich, 2013, pp. 760–765.
  * [14] Y. Mo and R. M. Murray, “Privacy preserving average consensus,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 2, pp. 753–765, 2017\.
  * [15] D. B. West, _Introduction to Graph Theory_. Upper Saddle River, New Jersey, 2001.
  * [16] A. I. Rikos, T. Charalambous, K. H. Johansson, and C. N. Hadjicostis, “Privacy-preserving event-triggered quantized average consensus,” in _Proceedings of $59^{th}$ IEEE Conference on Decision and Control_, 2020, pp. 6246–6253.
  * [17] ——, “Distributed event-triggered algorithms for finite-time privacy-preserving quantized average consensus,” _arXiv preprint arXiv:2102.06778_ , 2021.
  * [18] Q. Li, I. Cascudo, and M. G. Christensen, “Privacy-preserving distributed average consensus based on additive secret sharing,” in _European Signal Processing Conferenec (EUSIPCO)_ , 2019, pp. 1–5.
  * [19] C. N. Hadjicostis and A. D. Dominguez-Garcia, “Privacy-preserving distributed averaging via homomorphically encrypted ratio consensus,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 9, pp. 3887–3894, 2020.
  * [20] M. Ruan, H. Gao, and Y. Wang, “Secure and privacy-preserving consensus,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 10, pp. 4035–4049, 2019.
  * [21] C. N. Hadjicostis, “Privary preserving distributed average consensus via homomorphic encryption,” in _Proceedings of $57^{th}$ IEEE Conference on Decision and Control (CDC)_, 2018, pp. 1258–1263.
  * [22] T. Charalambous, N. E. Manitara, and C. N. Hadjicostis, “Privacy-preserving average consensus over digraphs in the presence of time delays,” in _Proceedings of $57^{th}$ Annual Allerton Conference on Communication, Control, and Computing_, 2019, pp. 238–245.
  * [23] L. Xiao and S. Boyd, “Fast linear iterations for distributed averaging,” _Systems and Control Letters_ , vol. 53, no. 1, pp. 65–78, September 2004\.
  * [24] N. E. Manitara and C. N. Hadjicostis, “Distributed stopping in average consensus via event-triggered strategies,” in _Proceedings of $51^{st}$ Annual Allerton Conference on Communication, Control, and Computing_, 2013, pp. 1336–1343.


# Peak height distributions of certain nonisotropic Gaussian random fields

Dan Cheng
Arizona State University

###### Abstract

The exact and explicit peak height distribution for smooth nonstationary unit-
variance Gaussian processes is derived such that it depends only on a
conditional correlation parameter. Further, we study the peak height
distributions of certain stationary nonisotropic Gaussian fields, including
planar Gaussian fields, cosine fields and anisotropic Gaussian fields, and
derive their exact formulas respectively. These extend the results in [6, 7]
from isotropy to nonisotropic cases.

Keywords: Gaussian random fields; Peak; Height distribution; Nonstationary;
Stationary; Isotropic; Nonisotropic; Anisotropic; Cosine fields; Planar
Gaussian fields.

Mathematics Subject Classification: 60G15, 60G60, 62G32, 15B52.

##  1 Introduction

The study of peaks (local maxima) of Gaussian random fields, especially of
their height distribution and expected number, is important in probability
theory [1, 9] and has applications in many areas such as statistics [8, 5, 4],
neuroimaging [12, 11] and astronomy [3]. However, the exact formulae for the
peak height distributions of Gaussian fields are difficult of obtain, mainly
due to the hard computation involving Hessian matrices. The only exception is
the isotropic Gaussian fields, where one can apply random matrices techniques
(GOE and GOI matrices) to deal with the Hessian and hence obtain the formalae
[6, 7]. In this paper, we investigate certain nonisotropic Gaussian fields to
obtain the exact peak height distribution formulae for a wider class of
Gaussian fields.

Specifically, let $X=\\{X(t),t\in T\\}$ be a centered, unit-variance, smooth
Gaussian random field living on a parameter set $T\subset{\mathbb{R}}^{N}$.
Throughout this paper, we assume that the filed $X$ satisfies the smoothness
condition in [6], which is slightly stronger than $C^{2}$ but can be implied
by $C^{3}$, and the regularity condition therein, such that the peak height
distribution can be computed by the Kac-Rice formula. We now introduce the
peak height distribution of $X$. It is convenient to consider the tail
probability, that is, the probability that the peak height exceeds a fixed
threshold at that point, conditioned on the event that the point is a peak
(local maximum) of $X$. Such conditional probability can be defined as

$F_{t}(u)={\mathbb{P}}[X(t)>u\ |\ t\text{ is a local maximum of }X].$

It has been shown in [6] that the peak height distribution of $X$ at $t$ is
given by

$F_{t}(u)=\frac{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{X(t)>u\\}}\mathbbm{1}_{\\{\nabla^{2}X(t)\prec
0\\}}|\nabla X(t)=0\right]}{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{\nabla^{2}X(t)\prec 0\\}}|\nabla
X(t)=0\right]},\quad u\in{\mathbb{R}},$ (1.1)

where $\nabla X(t)$ and $\nabla^{2}X(t)$ are respectively the gradient and
Hessian of $X$, and $\nabla^{2}X(t)\prec 0$ means that the Hessian is negative
definite. Taking the derivative, we obtain the peak height density and denote
it by $h_{t}(x)=-F_{t}(x)$, $x\in{\mathbb{R}}$.

In particular, if the field $X$ is stationary, then we simply write
$F_{t}(u)=F(u)$ and $h_{t}(x)=h(x)$ since they do not depend on the location
$t$. Moreover, due to the stationarity, $\nabla X(t)$ is independent of both
$X(t)$ and $\nabla^{2}X(t)$ for each fixed $t$, and hence the peak height
distribution in (1.1) becomes

$\begin{split}F(u)&=\frac{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{X(t)>u\\}}\mathbbm{1}_{\\{\nabla^{2}X(t)\prec
0\\}}\right]}{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{\nabla^{2}X(t)\prec 0\\}}\right]},\quad
u\in{\mathbb{R}}.\end{split}$ (1.2)

It can be seen from (1.1) or (1.2) that the expectation involving the
determinant and indicator function on the Hessian would be hard. We
investigate in this paper certain nonisotropic Gaussian fields with nice
properties that can simplify such expectation to derive the exact peak height
distributions. More specifically, the sections below will be focused on
nonstationary Gaussian processes, planar Gaussian fields, cosine fields and
anisotropic Gaussian fields in sequence.

Throughout this paper, we denote by $\phi(x)$ and $\Phi(x)$ the pdf and cdf of
the standard normal distribution $N(0,1)$, respectively. Let
$\Psi(x)=1-\Phi(x)$ be tail probability of $N(0,1)$. For a standard bivariate
normal vector $(Z_{1},Z_{2})$ with ${\rm Var}(Z_{1})={\rm Var}(Z_{2})=1$ and
${\rm Corr}(Z_{1},Z_{2})=\rho$, denote the joint density by

$\begin{split}P_{\rho}(x,y)=\frac{1}{2\pi\sqrt{1-\rho^{2}}}e^{-\frac{x^{2}-2\rho
xy+y^{2}}{2(1-\rho^{2})}},\end{split}$ (1.3)

and denote the expectation taken under the this density by

$\begin{split}{\mathbb{E}}_{\rho}[g(Z_{1},Z_{2})]=\int_{{\mathbb{R}}^{2}}g(z_{1},z_{2})P_{\rho}(z_{1},z_{2})dz_{1}dz_{2},\end{split}$
(1.4)

where $g$ is a measurable function. Denote by $I_{n}$ the identity matrix of
size $n$. For each $t\in{\mathbb{R}}^{N}$, write $t=(t_{1},\ldots,t_{N})$. For
$t,s\in{\mathbb{R}}^{N}$, denote their Euclidean distance by $\|t-s\|$.

##  2 Nonstationary Gaussian Processes

In this section, we study the case that $\\{X(t),t\in{\mathbb{R}}\\}$ is a
smooth nonstationary Gaussian process with variance 1. The main result is
stated as follows, showing that the peak height distribution depends only on a
single parameter $\rho(t)$.

######  Theorem 2.1.

Let $\\{X(t),t\in{\mathbb{R}}\\}$ be a smooth, unit-variance, nonstationary
Gaussian process. Then the peak height distribution and density are given
respectively by

$\begin{split}F_{t}(u)&=\Psi\left(\frac{u}{\sqrt{1-\rho^{2}(t)}}\right)-\sqrt{2\pi}\rho(t)\phi(u)\Phi\left(\frac{-\rho(t)u}{\sqrt{1-\rho^{2}(t)}}\right),\quad
u\in{\mathbb{R}},\\\
h_{t}(x)&=\sqrt{1-\rho^{2}(t)}\phi\left(\frac{x}{\sqrt{1-\rho^{2}(t)}}\right)-\sqrt{2\pi}\rho(t)x\phi(x)\Phi\left(\frac{-\rho(t)x}{\sqrt{1-\rho^{2}(t)}}\right),\quad
x\in{\mathbb{R}}.\end{split}$ (2.1)

Here $\rho(t)$ is the conditional correlation of $X(t)$ and
$X^{\prime\prime}(t)$ given $X^{\prime}(t)=0$,

$\begin{split}\rho(t)={\rm
Corr}[X(t),X^{\prime\prime}(t)|X^{\prime}(t)=0]=\frac{-\lambda_{1}(t)}{\sqrt{\lambda_{2}(t)-\frac{r^{2}(t)}{\lambda_{1}(t)}}},\end{split}$
(2.2)

where $\lambda_{1}(t)={\rm Var}(X^{\prime}(t))$, $\lambda_{2}(t)={\rm
Var}(X^{\prime\prime}(t))$ and
$r(t)={\mathbb{E}}[X^{\prime}(t)X^{\prime\prime}(t)]$.

###### Proof.

Note that ${\rm Var}(X(t))\equiv 1$, taking the first and second derivatives
on both sides yields ${\mathbb{E}}[X(t)X^{\prime}(t)]=0$ and ${\rm
Var}(X^{\prime}(t))=-{\mathbb{E}}[X(t)X^{\prime\prime}(t)]$. Then the random
vector $(X(t),X^{\prime}(t),X^{\prime\prime}(t))$ has the variance-covariance
matrix

$\begin{split}\begin{pmatrix}{\rm
Var}(X(t))&{\mathbb{E}}[X(t)X^{\prime}(t)]&{\mathbb{E}}[X(t)X^{\prime\prime}(t)]\\\
{\mathbb{E}}[X(t)X^{\prime}(t)]&{\rm
Var}(X^{\prime}(t))&{\mathbb{E}}[X^{\prime}(t)X^{\prime\prime}(t)]\\\
{\mathbb{E}}[X(t)X^{\prime\prime}(t)]&{\mathbb{E}}[X^{\prime}(t)X^{\prime\prime}(t)]&{\rm
Var}(X^{\prime\prime}(t))\end{pmatrix}=\begin{pmatrix}1&0&-\lambda_{1}(t)\\\
0&\lambda_{1}(t)&r(t)\\\
-\lambda_{1}(t)&r(t)&\lambda_{2}(t)\end{pmatrix}\end{split}.$

Since $X$ is a 1D Gaussian process, the peak height distribution in (1.1)
becomes

$\begin{split}F_{t}(u)&=\frac{{\mathbb{E}}\left[X^{\prime\prime}(t)\mathbbm{1}_{\\{X^{\prime\prime}(t)<0\\}}\mathbbm{1}_{\\{X(t)>u\\}}|X^{\prime}(t)=0\right]}{{\mathbb{E}}\left[X^{\prime\prime}(t)\mathbbm{1}_{\\{X^{\prime\prime}(t)<0\\}}|X^{\prime}(t)=0\right]}.\end{split}$
(2.3)

Note that

$\begin{split}{\rm Var}[X(t)|X^{\prime}(t)=0]&=1,\\\ {\rm
Var}[X^{\prime\prime}(t)|X^{\prime}(t)=0]&=\lambda_{2}(t)-\frac{r^{2}(t)}{\lambda_{1}(t)}:=\delta^{2}(t),\\\
{\rm
Cov}[X(t),X^{\prime\prime}(t)|X^{\prime}(t)=0]&=-\lambda_{1}(t),\end{split}$

implying

$\begin{split}(X(t),X^{\prime\prime}(t)|X^{\prime}(t)=0)\sim
N\left(\begin{pmatrix}0\\\
0\end{pmatrix},\,\begin{pmatrix}1&-\lambda_{1}(t)\\\
-\lambda_{1}(t)&\delta^{2}(t)\end{pmatrix}\right).\end{split}$ (2.4)

In particular, the correlation of $X(t)$ and $X^{\prime\prime}(t)$ conditional
on $X^{\prime}(t)=0$ is given by (2.2). Since
$(X^{\prime\prime}(t)|X^{\prime}(t)=0)\sim N(0,\,\delta^{2}(t))$, the
denominator in (2.3) becomes

$\begin{split}{\mathbb{E}}\left[X^{\prime\prime}(t)\mathbbm{1}_{\\{X^{\prime\prime}(t)<0\\}}|X^{\prime}(t)=0\right]=\delta(t)\int_{-\infty}^{0}y\phi(y)dy=\frac{-\delta(t)}{\sqrt{2\pi}}.\end{split}$
(2.5)

On the other hand, applying (2.4) and (2.2), as well as integration by parts,
the numerator in (2.3) becomes

$\begin{split}&\quad{\mathbb{E}}\left[X^{\prime\prime}(t)\mathbbm{1}_{\\{X^{\prime\prime}(t)<0\\}}\mathbbm{1}_{\\{X(t)>u\\}}|X^{\prime}(t)=0\right]\\\
&=\delta(t)\int_{u}^{\infty}dz\int_{-\infty}^{0}\frac{w}{2\pi\sqrt{1-\rho^{2}(t)}}e^{-\frac{z^{2}-2\rho(t)zw+w^{2}}{2(1-\rho^{2}(t))}}dw\\\
&=\delta(t)(1-\rho^{2}(t))\int_{\frac{u}{\sqrt{1-\rho^{2}(t)}}}^{\infty}dx\int_{-\infty}^{0}\frac{y}{2\pi}e^{-\frac{x^{2}-2\rho(t)xy+y^{2}}{2}}dy\\\
&=\delta(t)(1-\rho^{2}(t))\int_{\frac{u}{\sqrt{1-\rho^{2}(t)}}}^{\infty}\frac{1}{\sqrt{2\pi}}e^{-\frac{(1-\rho^{2}(t))x^{2}}{2}}dx\int_{-\infty}^{0}\frac{y}{\sqrt{2\pi}}e^{-\frac{(y-\rho(t)x)^{2}}{2}}dy\\\
&=\delta(t)(1-\rho^{2}(t))\int_{\frac{u}{\sqrt{1-\rho^{2}(t)}}}^{\infty}\frac{1}{\sqrt{2\pi}}e^{-\frac{(1-\rho^{2}(t))x^{2}}{2}}[\rho(t)x\Phi(-\rho(t)x)-\phi(\rho(t)x)]dx\\\
&=-\frac{\delta(t)}{\sqrt{2\pi}}\left(\int_{\frac{u}{\sqrt{1-\rho^{2}(t)}}}^{\infty}\rho(t)\Phi(-\rho(t)x)de^{-\frac{(1-\rho^{2}(t))x^{2}}{2}}+(1-\rho^{2}(t))\Psi\left(\frac{u}{\sqrt{1-\rho^{2}(t)}}\right)\right)\\\
&=-\frac{\delta(t)}{\sqrt{2\pi}}\left(\Psi\left(\frac{u}{\sqrt{1-\rho^{2}(t)}}\right)-\sqrt{2\pi}\rho(t)\phi(u)\Phi\left(-\frac{\rho(t)u}{\sqrt{1-\rho^{2}(t)}}\right)\right),\end{split}$
(2.6)

where we have used change of variables $x=\frac{z}{\sqrt{1-\rho^{2}(t)}}$ and
$y=\frac{w}{\sqrt{1-\rho^{2}(t)}}$ in the third line. Taking the ratio between
(2.6) and (2.5), we obtain the peak height distribution

$F_{t}(u)=\Psi\left(\frac{u}{\sqrt{1-\rho^{2}(t)}}\right)-\sqrt{2\pi}\rho(t)\phi(u)\Phi\left(-\frac{\rho(t)u}{\sqrt{1-\rho^{2}(t)}}\right),\quad
u\in{\mathbb{R}}.$

Further, taking the derivative, we obtain the peak height density

$\begin{split}h_{t}(x)&=\frac{1}{\sqrt{1-\rho^{2}(t)}}\phi\left(\frac{x}{\sqrt{1-\rho^{2}(t)}}\right)-\frac{\sqrt{2\pi}\rho^{2}(t)}{\sqrt{1-\rho^{2}(t)}}\phi(x)\phi\left(-\frac{\rho(t)x}{\sqrt{1-\rho^{2}(t)}}\right)\\\
&\quad-\sqrt{2\pi}\rho(t)x\phi(x)\Phi\left(-\frac{\rho(t)x}{\sqrt{1-\rho^{2}(t)}}\right)\\\
&=\sqrt{1-\rho^{2}(t)}\phi\left(\frac{x}{\sqrt{1-\rho^{2}(t)}}\right)-\sqrt{2\pi}\rho(t)x\phi(x)\Phi\left(-\frac{\rho(t)x}{\sqrt{1-\rho^{2}(t)}}\right),\end{split}$

completing the proof. ∎

Here we make some remarks on the comparison between the peak height
distributions for nonstationary and stationary Gaussian processes. Let $Z$ be
a 1D stationary Gaussian process with covariance
${\mathbb{E}}[Z(t)Z(s)]=\varphi\left(|t-s|^{2}\right)$, where
$t,s\in{\mathbb{R}}$ and $\varphi:[0,\infty)\to{\mathbb{R}}$. The peak height
distribution and density can be found in [9] or be derived from the general
formula for Gaussian fields in [6, 7] such that

$\begin{split}F(u)&=\Psi\left(\frac{x}{\sqrt{1-\kappa^{2}/3}}\right)+\frac{\sqrt{2\pi}\kappa}{\sqrt{3}}\phi(x)\Phi\left(\frac{\kappa
x}{\sqrt{3-\kappa^{2}}}\right),\\\
h(x)&=\sqrt{1-\kappa^{2}/3}\phi\left(\frac{x}{\sqrt{1-\kappa^{2}/3}}\right)+\frac{\sqrt{2\pi}\kappa}{\sqrt{3}}x\phi(x)\Phi\left(\frac{\kappa
x}{\sqrt{3-\kappa^{2}}}\right),\end{split}$ (2.7)

where $\kappa=-\varphi^{\prime}(0)/\sqrt{\varphi^{\prime\prime}(0)}$, or
equivalently $\kappa=\sqrt{3}{\rm Var}[Z^{\prime}(t)]/\sqrt{{\rm
Var}[Z^{\prime\prime}(t)]}$ since ${\rm
Var}[Z^{\prime}(t)]=-2\varphi^{\prime}(0)$ and ${\rm
Var}[Z^{\prime\prime}(t)]=12\varphi^{\prime\prime}(0)$ as shown in [6]. Note
that $r(t)=0$ due to the independence of $X^{\prime}(t)$ and
$X^{\prime\prime}(t)$ under stationarity, thus $\kappa=-\sqrt{3}\rho(t)$ by
(2.2). Plugging this into (2.7), we see that the two peak height distributions
coincide for the stationary case.

As a nice and important property for the peak height distribution of
nonstationary Gaussian processes in (2.1), we see that the formula depends
only on a single parameter $\rho(t)$, which is the conditional correlation of
$X(t)$ and $X^{\prime\prime}(t)$ given $X^{\prime}(t)=0$. It has been shown in
[6, 7] that, for isotropic Gaussian fields in ${\mathbb{R}}^{N}$, the peak
height distribution also depends only on one parameter. However, this property
usually does not hold even for stationary Gaussian fields; see Section 3 below
on the case of 2D stationary Gaussian fields.

The peak height distribution in (2.7) has been used in computing p-values in
statistical applications such as peak/signal detection [8] and change point
detection [5]. We are expecting to apply the newly derived formula in (2.1) to
perform multiple testings for both signal and change point detection under
nonstationary Gaussian noise. Particularly, the nice property that the peak
height distribution depends only on a single parameter would be a big
advantage in parameter estimation and statistical inference.

##  3 Planar Stationary Gaussian Random Fields

In this section, we study stationary (nonisotropic) Gaussian fields on
${\mathbb{R}}^{2}$. Generally, the peak height distribution is hard to
evaluate due to the complexity of the gradient and Hessian matrix. We shall
add certain independence assumptions on the gradient and Hessian (see (3.2)
below) to obtain the peak height distribution formula in Theorem 3.2.

###  3.1 Notations and Assumptions

Let $\\{X(t):t\in{\mathbb{R}}^{2}\\}$ be a smooth stationary planar Gaussian
random field (i.e., parameterized on ${\mathbb{R}}^{2}$). Let

$X_{i}(t)=\frac{\partial X}{\partial t_{i}}(t)\quad\text{\rm and }\quad
X_{kl}(t)=\frac{\partial^{2}X}{\partial t_{k}\partial t_{l}}(t),\quad
i,k,l=1,2.$

Note that the distribution of $(X(t),X_{i}(t),X_{kl}(t))$ does not depend on
$t$ due to the stationarity of $X$. We have the following spectral
representations (see page 112 in [1]) for the covariance of the gradient and
Hessian,

${\mathbb{E}}[X_{i}(t)X_{j}(t)]=\int_{{\mathbb{R}}^{2}}x_{i}x_{j}\nu(dx)\quad\text{\rm
and
}\quad{\mathbb{E}}[X_{ij}(t)X_{kl}(t)]=\int_{{\mathbb{R}}^{4}}x_{i}x_{j}x_{k}x_{l}\nu(dx),$
(3.1)

where $\nu(\cdot)$ is a spectral measure. As mentioned before, for general
stationary Gaussian fields, it would be hard to compute the peak height
distribution. Here, in this section, we assume

${\mathbb{E}}[X_{1}(t)X_{2}(t)]={\mathbb{E}}[X_{11}(t)X_{12}(t)]={\mathbb{E}}[X_{22}(t)X_{12}(t)]=0.$
(3.2)

The following example of covariance $C(t)={\mathbb{E}}[X(t)X(0)]$, where
$t=(t_{1},t_{2})\in{\mathbb{R}}^{2}$, is of this case:

$C(t_{1},t_{2})=C(-t_{1},t_{2})\quad\text{ and }\quad
C(t_{1},t_{2})=C(t_{1},-t_{2}).$

Equivalently, the spectral measure $\nu$ in (3.1) is symmetric in each
direction $\lambda_{1}$ and $\lambda_{2}$, separately. This is called quadrant
symmetry. Recall that an isotropic Gaussian field $Y$ has covariance such that
${\mathbb{E}}[Y(t)Y(s)]$ depends only on $\|t-s\|$. Thus such quadrant
symmetry property is weaker than isotropy.

Here, ${\mathbb{E}}[X_{1}(t)X_{2}(t)]=0$ indicates that ${\rm Cov}[\nabla
X(t)]$ is a diagonal matrix. Note that, by (3.1),
${\mathbb{E}}[X_{11}(t)X_{22}(t)]={\rm Var}[X_{12}(t)]$. Hence, under our
assumption (3.2), the covariance of the Hessian can be written as

$\begin{split}{\rm
Cov}[X_{11}(t),X_{22}(t),X_{12}(t)]=\begin{pmatrix}\sigma_{1}^{2}&\sigma_{3}^{2}&0\\\
\sigma_{3}^{2}&\sigma_{2}^{2}&0\\\ 0&0&\sigma_{3}^{2}\\\
\end{pmatrix}.\end{split}$

###  3.2 Peak Height Distribution of Planar Stationary Gaussian Fields

We first derive the following result that would be helpful in simplifying the
peak height distribution.

######  Lemma 3.1.

Let $\\{X(t),t\in{\mathbb{R}}^{2}\\}$ be a smooth, unit-variance, stationary
Gaussian process satisfying ${\rm Cov}[X_{1}(t),X_{2}(t)]={\rm
diag}(\gamma_{1}^{2},\gamma_{2}^{2})$ and
${\mathbb{E}}[X_{11}(t)X_{12}(t)]={\mathbb{E}}[X_{22}(t)X_{12}(t)]=0$, where
$\gamma_{1},\gamma_{2}>0$. Let
$Z(t_{1},t_{2})=X(t_{1}/\gamma_{1},t_{2}/\gamma_{2})$. Then $Z$ is stationary
with ${\rm Cov}[Z_{1}(t),Z_{2}(t)]=I_{2}$ and
${\mathbb{E}}[Z_{11}(t)Z_{12}(t)]={\mathbb{E}}[Z_{22}(t)Z_{12}(t)]=0$, and it
has the same peak height distribution as that of $X$.

###### Proof.

First, note that the transformation is linear, and thus keeps the stationarity
of $Z$. Since $Z(t_{1},t_{2})=X(t_{1}/\gamma_{1},t_{2}/\gamma_{2})$, one has

$\begin{split}Z_{i}(t_{1},t_{2})&=X_{i}(t_{1}/\gamma_{1},t_{2}/\gamma_{2})/\gamma_{i},\quad
i=1,\\\
Z_{jk}(t_{1},t_{2})&=X_{jk}(t_{1}/\gamma_{1},t_{2}/\gamma_{2})/(\gamma_{j}\gamma_{k}),\quad
j,k=1,2.\end{split}$

Therefore,

$\begin{split}{\rm Cov}[Z_{1}(t),Z_{2}(t)]&={\rm
Cov}[X_{1}(t)/\gamma_{1},X_{2}(t)/\gamma_{2}]=I_{2},\\\
{\mathbb{E}}[Z_{ii}(t)Z_{ij}(t)]&={\mathbb{E}}[X_{ii}(t)X_{ij}(t)]/\left(\gamma_{i}^{3}\gamma_{j}\right)=0,\quad
i\neq j.\end{split}$

Hence the peak height distribution in (1.2) becomes

$\begin{split}F(u)&=\frac{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}Z(t_{1},t_{2}))|\mathbbm{1}_{\\{Z(t_{1},t_{2})>u\\}}\mathbbm{1}_{\\{\nabla^{2}Z(t_{1},t_{2})\prec
0\\}}\right]}{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}Z(t_{1},t_{2}))|\mathbbm{1}_{\\{\nabla^{2}Z(t_{1},t_{2})\prec
0\\}}\right]}\\\ &=\frac{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t_{1}/\gamma_{1},t_{2}/\gamma_{2}))|\mathbbm{1}_{\\{X(t_{1}/\gamma_{1},t_{2}/\gamma_{2})>u\\}}\mathbbm{1}_{\\{A\nabla^{2}X(t_{1}/\gamma_{1},t_{2}/\gamma_{2})A\prec
0\\}}\right]/(\gamma_{1}^{2}\gamma_{2}^{2})}{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t_{1}/\gamma_{1},t_{2}/\gamma_{2}))|\mathbbm{1}_{\\{A\nabla^{2}X(t_{1}/\gamma_{1},t_{2}/\gamma_{2})A\prec
0\\}}\right]/(\gamma_{1}^{2}\gamma_{2}^{2})}\\\
&=\frac{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t_{1}/\gamma_{1},t_{2}/\gamma_{2}))|\mathbbm{1}_{\\{X(t_{1}/\gamma_{1},t_{2}/\gamma_{2})>u\\}}\mathbbm{1}_{\\{\nabla^{2}X(t_{1}/\gamma_{1},t_{2}/\gamma_{2})\prec
0\\}}\right]}{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t_{1}/\gamma_{1},t_{2}/\gamma_{2}))|\mathbbm{1}_{\\{\nabla^{2}X(t_{1}/\gamma_{1},t_{2}/\gamma_{2})\prec
0\\}}\right]},\end{split}$

where $A={\rm diag}(1/\gamma_{1},1/\gamma_{2})$. This implies that $Z$ and $X$
have the same peak height distribution. ∎

By Lemma 3.1, for the general case with ${\rm Cov}[X_{1},X_{2}]={\rm
diag}(\gamma_{1}^{2},\gamma_{2}^{2})$, we can study instead
$Z(t_{1},t_{2})=X(t_{1}/\gamma_{1},t_{2}/\gamma_{2})$ which satisfies ${\rm
Cov}[Z_{1},Z_{2}]=I_{2}$. Thus, without loss of generality, let us assume
${\rm Cov}[X_{1},X_{2}]=I_{2}$. Define the function

$\begin{split}&\quad g(z_{1},z_{2},x|a_{1}^{2},a_{2}^{2},a_{3}^{2})\\\
&=\alpha_{1}a_{2}\Bigg{\\{}\left[\left(z_{1}-\frac{x}{a_{1}}\right)\left(z_{2}-\frac{x}{a_{2}}\right)-b\right]\left[\Phi\left(\sqrt{\left(z_{1}-\frac{x}{a_{1}}\right)\left(z_{2}-\frac{x}{a_{2}}\right)/b}\right)-\frac{1}{2}\right]\\\
&\quad+\sqrt{b\left(z_{1}-\frac{x}{a_{1}}\right)\left(z_{2}-\frac{x}{a_{2}}\right)}\phi\left(\sqrt{\left(z_{1}-\frac{x}{a_{1}}\right)\left(z_{2}-\frac{x}{a_{2}}\right)/b}\right)\Bigg{\\}}\mathbbm{1}_{\left\\{z_{1}<\frac{x}{a_{1}}\right\\}}\mathbbm{1}_{\left\\{z_{2}<\frac{x}{a_{2}}\right\\}},\end{split}$
(3.3)

where $a_{1},a_{2},a_{3}>0$ and $b=a_{3}^{2}/(a_{1}a_{2})$. The peak height
distribution for 2D stationary Gaussian fields is derived as follows.

######  Theorem 3.2.

Let $\\{X(t),t\in{\mathbb{R}}^{2}\\}$ be a smooth, unit-variance, stationary
Gaussian process with ${\rm Cov}[X_{1},X_{2}]=I_{2}$ and
${\mathbb{E}}[X_{11}X_{12}]={\mathbb{E}}[X_{22}X_{12}]=0$. Let ${\rm
Var}(X_{11})=\sigma_{1}^{2}$, ${\rm Var}(X_{22})=\sigma_{2}^{2}$ and ${\rm
Var}(X_{12})=\sigma_{3}^{2}$. Then the peak height density is given by

$\begin{split}h(x)&=\frac{\phi(x){\mathbb{E}}_{\tilde{\rho}}[g(\tilde{Z}_{1},\tilde{Z}_{2},x|\sigma_{1}^{2}-1,\sigma_{2}^{2}-1,\sigma_{3}^{2})]}{{\mathbb{E}}_{\rho}[g(Z_{1},Z_{2},0|\sigma_{1}^{2},\sigma_{2}^{2},\sigma_{3}^{2})]},\end{split}$
(3.4)

where $(Z_{1},Z_{2})$ and $(\tilde{Z}_{1},\tilde{Z}_{2})$ are respectively
standard bivariate normal vectors with correlations

$\begin{split}\rho=\frac{\sigma_{3}^{2}}{\sigma_{1}\sigma_{2}}\quad\text{\rm
and
}\quad\tilde{\rho}=\frac{\sigma_{3}^{2}-1}{\sqrt{(\sigma_{1}^{2}-1)(\sigma_{1}^{2}-1)}},\end{split}$

${\mathbb{E}}_{\rho}$ and ${\mathbb{E}}_{\tilde{\rho}}$ are the corresponding
bivariate normal expectations defined in (1.4), and the function $g$ is
defined in (3.3).

###### Proof.

Note that, the peak height distribution in (1.2) for stationary Gaussian
fields can be written as

$\begin{split}F(u)&=\frac{\int_{u}^{\infty}\phi(x){\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{\nabla^{2}X(t)\prec
0\\}}|X(t)=x\right]dx}{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{\nabla^{2}X(t)\prec
0\\}}\right]}.\end{split}$

Taking the derivative, we obtain the peak height density

$\begin{split}h(x)=\frac{\phi(x){\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{\nabla^{2}X(t)\prec
0\\}}|X(t)=x\right]}{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{\nabla^{2}X(t)\prec
0\\}}\right]}.\end{split}$ (3.5)

It is well known that, for a symmetric $2\times 2$ matrix, it is negative
definite if and only if both diagonal entries are negative and the determinant
is positive. Using this fact, we can write the denominator in (3.5) as

$\begin{split}J_{1}&={\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{\nabla^{2}X(t)\prec 0\\}}\right]\\\
&={\mathbb{E}}\left[(X_{11}(t)X_{22}(t)-X_{12}^{2}(t))\mathbbm{1}_{\\{X_{11}(t)<0\\}}\mathbbm{1}_{\\{X_{22}(t)<0\\}}\mathbbm{1}_{\\{X_{11}(t)X_{22}(t)-X_{12}^{2}(t)>0\\}}\right].\end{split}$

Let $Z_{1}=X_{11}(t)/\sigma_{1}$, $Z_{2}=X_{22}(t)/\sigma_{2}$ and
$Z_{3}=X_{12}(t)/\sigma_{3}$. Note that, by the assumption, $Z_{3}$ is
independent of both $Z_{1}$ and $Z_{2}$, thus we have

$\begin{split}J_{1}&={\mathbb{E}}\left[(\sigma_{1}\sigma_{2}Z_{1}Z_{2}-\sigma_{3}^{2}Z_{3}^{2})\mathbbm{1}_{\\{Z_{1}<0\\}}\mathbbm{1}_{\\{Z_{2}<0\\}}\mathbbm{1}_{\\{\sigma_{1}\sigma_{2}Z_{1}Z_{2}-\sigma_{3}^{2}Z_{3}^{2}>0\\}}\right].\\\
&=\int_{-\infty}^{0}\int_{-\infty}^{0}(\sigma_{1}\sigma_{2}z_{1}z_{2}-\sigma_{3}^{2}z_{3}^{2})\left(\int_{-\sqrt{\frac{z_{1}z_{2}}{\rho}}}^{\sqrt{\frac{z_{1}z_{2}}{\rho}}}\phi(z_{3})dz_{3}\right)P_{\rho}(z_{1},z_{2})dz_{1}dz_{2}\\\
&=2\sigma_{1}\sigma_{2}\int_{-\infty}^{0}\int_{-\infty}^{0}(z_{1}z_{2}-\rho
z_{3}^{2})\left(\int_{0}^{\sqrt{\frac{z_{1}z_{2}}{\rho}}}\phi(z_{3})dz_{3}\right)P_{\rho}(z_{1},z_{2})dz_{1}dz_{2},\end{split}$

where $\rho=\sigma_{3}^{2}/(\sigma_{1}\sigma_{2})$ and $P_{\rho}(z_{1},z_{2})$
is the joint density of a standard bivariate normal vector with correlation
$\rho$ defined in (1.3). Note that

$\int_{0}^{u}\phi(z_{3})dz_{3}=\Phi(u)-\frac{1}{2},\quad\int_{0}^{u}z_{3}^{2}\phi(z_{3})dz_{3}=-u\phi(u)+\Phi(u)-\frac{1}{2},$
(3.6)

we obtain

$\begin{split}J_{1}&=2\sigma_{1}\sigma_{2}\int_{-\infty}^{0}\int_{-\infty}^{0}\left[(z_{1}z_{2}-\rho)\left(\Phi\left(\sqrt{\frac{z_{1}z_{2}}{\rho}}\right)-\frac{1}{2}\right)+\sqrt{\rho
z_{1}z_{2}}\phi\left(\sqrt{\frac{z_{1}z_{2}}{\rho}}\right)\right]P_{\rho}(z_{1},z_{2})dz_{1}dz_{2}\\\
&=2{\mathbb{E}}_{\rho}[g(Z_{1},Z_{2},0|\sigma_{1}^{2},\sigma_{2}^{2},\sigma_{3}^{2})],\end{split}$
(3.7)

where the function $g$ is defined in (3.3).

Now, we turn to computing the numerator in (3.5). Note that
${\mathbb{E}}[X(t)\nabla^{2}X(t)]=-{\rm Cov}[\nabla X(t)]=-I_{2}$. Applying
the conditional formula for multivariate normal variables, we have

$\begin{split}{\mathbb{E}}[\nabla^{2}X(t)|X(t)=x]&=-I_{2}x,\\\ {\rm
Cov}[X_{11}(t),X_{22}(t),X_{12}(t)|X(t)=x]&=\begin{pmatrix}\sigma_{1}^{2}-1&\sigma_{3}^{2}-1&0\\\
\sigma_{3}^{2}-1&\sigma_{2}^{2}-1&0\\\ 0&0&\sigma_{3}^{2}\\\
\end{pmatrix},\end{split}$

implying that

$\begin{split}(X_{11}(t),X_{22}(t),X_{12}(t)|X(t)=x)\sim
N\left(\begin{pmatrix}-x\\\ -x\\\ 0\\\
\end{pmatrix},\begin{pmatrix}\sigma_{1}^{2}-1&\sigma_{3}^{2}-1&0\\\
\sigma_{3}^{2}-1&\sigma_{2}^{2}-1&0\\\ 0&0&\sigma_{3}^{2}\\\
\end{pmatrix}\right).\end{split}$

Let $\tilde{Z}_{1}=(X_{11}(t)+x)/(\sqrt{\sigma_{1}^{2}-1})$,
$\tilde{Z}_{2}=(X_{22}(t)+x)/(\sqrt{\sigma_{2}^{2}-1})$ and
$\tilde{Z}_{3}=X_{12}(t)/\sigma_{3}$. Then

$\begin{split}(\tilde{Z}_{1},\tilde{Z}_{2},\tilde{Z}_{3}|X(t)=x)\sim
N\left(\begin{pmatrix}0\\\ 0\\\ 0\\\
\end{pmatrix},\begin{pmatrix}1&\tilde{\rho}&0\\\ \tilde{\rho}&1&0\\\ 0&0&1\\\
\end{pmatrix}\right),\end{split}$

where

$\tilde{\rho}=\frac{\sigma_{3}^{2}-1}{\sqrt{(\sigma_{1}^{2}-1)(\sigma_{2}^{2}-1)}}.$

Now we have

$\begin{split}J_{2}&={\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{\nabla^{2}X(t)\prec
0\\}}|X(t)=x\right]\\\
&={\mathbb{E}}\left[(X_{11}(t)X_{22}(t)-X_{12}^{2}(t))\mathbbm{1}_{\\{X_{11}(t)<0\\}}\mathbbm{1}_{\\{X_{22}(t)<0\\}}\mathbbm{1}_{\\{X_{11}(t)X_{22}(t)-X_{12}^{2}(t)>0\\}}\big{|}X(t)=x\right]\\\
&={\mathbb{E}}\bigg{[}\left(\left(\sqrt{\sigma_{1}^{2}-1}\tilde{Z}_{1}-x\right)\left(\sqrt{\sigma_{2}^{2}-1}\tilde{Z}_{2}-x\right)-\sigma_{3}^{2}\tilde{Z}_{3}^{2}\right)\mathbbm{1}_{\left\\{\sqrt{\sigma_{1}^{2}-1}\tilde{Z}_{1}<x\right\\}}\mathbbm{1}_{\left\\{\sqrt{\sigma_{2}^{2}-1}\tilde{Z}_{2}<x\right\\}}\\\
&\qquad\times\mathbbm{1}_{\left\\{\left(\sqrt{\sigma_{1}^{2}-1}\tilde{Z}_{1}-x\right)\left(\sqrt{\sigma_{2}^{2}-1}\tilde{Z}_{2}-x\right)-\sigma_{3}^{2}\tilde{Z}_{3}^{2}>0\right\\}}\Big{|}X(t)=x\bigg{]}\\\
&=2\int_{-\infty}^{\frac{x}{\sqrt{\sigma_{1}^{2}-1}}}\int_{-\infty}^{\frac{x}{\sqrt{\sigma_{2}^{2}-1}}}\left(\left(\sqrt{\sigma_{1}^{2}-1}z_{1}-x\right)\left(\sqrt{\sigma_{2}^{2}-1}z_{2}-x\right)-\sigma_{3}^{2}z_{3}^{2}\right)\\\
&\qquad\times\int_{0}^{\frac{\sqrt{\left(\sqrt{\sigma_{1}^{2}-1}z_{1}-x\right)\left(\sqrt{\sigma_{2}^{2}-1}z_{2}-x\right)}}{\sigma_{3}}}\phi(z_{3})dz_{3}P_{\tilde{\rho}}(z_{1},z_{2})dz_{1}dz_{2}\\\
&=2\sqrt{\sigma_{1}^{2}-1}\sqrt{\sigma_{2}^{2}-1}\int_{-\infty}^{\frac{x}{\sqrt{\sigma_{1}^{2}-1}}}\int_{-\infty}^{\frac{x}{\sqrt{\sigma_{2}^{2}-1}}}\left(\left(z_{1}-\frac{x}{\sqrt{\sigma_{1}^{2}-1}}\right)\left(z_{2}-\frac{x}{\sqrt{\sigma_{2}^{2}-1}}\right)-bz_{3}^{2}\right)\\\
&\qquad\times\int_{0}^{\sqrt{\left(z_{1}-\frac{x}{\sqrt{\sigma_{1}^{2}-1}}\right)\left(z_{2}-\frac{x}{\sqrt{\sigma_{2}^{2}-1}}\right)/b}}\phi(z_{3})dz_{3}P_{\tilde{\rho}}(z_{1},z_{2})dz_{1}dz_{2},\end{split}$

where

$b=\frac{\sigma_{3}^{2}}{\sqrt{(\sigma_{1}^{2}-1)(\sigma_{2}^{2}-1)}}.$

It follows from (3.6) that we obtain

$\begin{split}J_{2}&=2\sqrt{\sigma_{1}^{2}-1}\sqrt{\sigma_{2}^{2}-1}\int_{-\infty}^{\frac{x}{\sqrt{\sigma_{1}^{2}-1}}}\int_{-\infty}^{\frac{x}{\sqrt{\sigma_{2}^{2}-1}}}\Bigg{[}\left(\left(z_{1}-\frac{x}{\sqrt{\sigma_{1}^{2}-1}}\right)\left(z_{2}-\frac{x}{\sqrt{\sigma_{2}^{2}-1}}\right)-b\right)\\\
&\qquad\times\left(\Phi\left(\sqrt{\left(z_{1}-\frac{x}{\sqrt{\sigma_{1}^{2}-1}}\right)\left(z_{2}-\frac{x}{\sqrt{\sigma_{2}^{2}-1}}\right)/b}\right)-\frac{1}{2}\right)\\\
&\qquad+\sqrt{b\left(z_{1}-\frac{x}{\sqrt{\sigma_{1}^{2}-1}}\right)\left(z_{2}-\frac{x}{\sqrt{\sigma_{2}^{2}-1}}\right)}\\\
&\qquad\times\phi\left(\sqrt{\left(z_{1}-\frac{x}{\sqrt{\sigma_{1}^{2}-1}}\right)\left(z_{2}-\frac{x}{\sqrt{\sigma_{2}^{2}-1}}\right)/b}\right)\Bigg{]}P_{\tilde{\rho}}(z_{1},z_{2})dz_{1}dz_{2}\\\
&=2{\mathbb{E}}_{\tilde{\rho}}[g(\tilde{Z}_{1},\tilde{Z}_{2},x|\sigma_{1}^{2}-1,\sigma_{2}^{2}-1,\sigma_{3}^{2})].\end{split}$
(3.8)

Plugging (3.7) and (3.8) into (3.5), we obtain the desired peak height density
in (3.4). ∎

##  4 Cosine Fields

The cosine random field on ${\mathbb{R}}^{N}$ is defined by (see [1] for
example)

$X(t)=\frac{1}{\sqrt{N}}\sum_{k=1}^{N}[\zeta_{k}\cos(\omega_{k}t_{k})+\zeta^{\prime}_{k}\sin(\omega_{k}t_{k})],$
(4.1)

where $\zeta_{k}$ and $\zeta^{\prime}_{k}$ are independent, standard Gaussian
variables and the $\omega_{k}$ are positive constants. It is a stationary,
nonisotropic, unit-variance, smooth Gaussian field with covariance

$C(t)={\mathbb{E}}[X(t)X(0)]=\frac{1}{N}\sum_{k=1}^{N}\cos(\omega_{k}t_{k}).$

Similarly to the proof of Lemma 3.1, the peak height distribution will not
change if we replace all $\omega_{k}$ by 1. Therefore, without loss of
generality, we assume $\omega_{k}=1$ for all $k=1,\ldots,N$ in the discussion
below.

######  Theorem 4.1.

Let $\\{X(t),t\in{\mathbb{R}}^{N}\\}$ be the cosine field defined in (4.1).
Then the peak height distribution is given by

$\begin{split}F(u)&=(2\pi)^{N/2}{\mathbb{E}}\left[\left(\prod_{i=1}^{N}|Z_{i}|\mathbbm{1}_{\\{Z_{i}<0\\}}\right)\mathbbm{1}_{\left\\{\sum_{i=1}^{N}Z_{i}\leq-\sqrt{N}u\right\\}}\right],\quad
u\geq 0,\end{split}$ (4.2)

where $Z_{1},\ldots,Z_{N}$ are independent standard normal variables.

###### Proof.

Let $\omega_{k}=1$ for all $k=1,\ldots,N$. We have the following covariance
structure

$\begin{split}{\rm Cov}(\nabla X(t))&=\frac{1}{N}I_{N},\\\ {\rm
Var}(X_{ii}(t))&=\frac{1}{N},\quad i=1,\ldots,N,\\\
{\mathbb{E}}[X_{ii}(t)X_{jj}(t)]&=0,\quad i\neq j.\end{split}$ (4.3)

Moreover,

$\begin{split}-\sum_{i=1}^{N}X_{ii}(t)=X(t)\quad\text{\rm and }\quad
X_{ij}(t)&=0,\quad i\neq j.\end{split}$ (4.4)

It follows from (4.3) and (4.4) that

$\begin{split}{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{\nabla^{2}X(t)\prec
0\\}}\right]&=(-1)^{N}{\mathbb{E}}\left[\prod_{i=1}^{N}X_{ii}(t)\mathbbm{1}_{\\{X_{ii}(t)<0\\}}\right]\\\
&=(-1)^{N}\left({\mathbb{E}}\left[X_{11}(t)\mathbbm{1}_{\\{X_{11}(t)<0\\}}\right]\right)^{N}=\frac{1}{(2\pi
N)^{N/2}}\end{split}$

and

$\begin{split}&\quad{\mathbb{E}}\left[|\text{det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{X(t)\geq
u\\}}\mathbbm{1}_{\\{\text{index}(\nabla^{2}X(t))=N\\}}\right]\\\
&=(-1)^{N}{\mathbb{E}}\left[\left(\prod_{i=1}^{N}X_{ii}(t)\mathbbm{1}_{\\{X_{ii}(t)<0\\}}\right)\mathbbm{1}_{\left\\{\sum_{i=1}^{N}X_{ii}(t)\leq-u\right\\}}\right]\\\
&=\frac{(-1)^{N}}{N^{N/2}}{\mathbb{E}}\left[\left(\prod_{i=1}^{N}Z_{i}\mathbbm{1}_{\\{Z_{i}<0\\}}\right)\mathbbm{1}_{\left\\{\sum_{i=1}^{N}Z_{i}\leq-\sqrt{N}u\right\\}}\right],\end{split}$

where $Z_{1},\ldots,Z_{N}$ are independent standard normal variables. The
desired result then follows from taking the ratio between the two expectations
above. ∎

##  5 Anisotropic Gaussian Random Fields

In this section, we study anisotropic Gaussian fields, which are useful models
in spatial statistics [2] and medical imaging [13].

###  5.1 Gaussian Orthogonally Invariant (GOI) Matrices

It is introduced in [7] that, an $N\times N$ random matrix $M=(M_{ij})_{1\leq
i,j\leq N}$ is called _Gaussian Orthogonally Invariant_ (GOI) with _covariance
parameter_ $c$, denoted by ${\rm GOI}(c)$, if it is symmetric and all entries
are centered Gaussian variables such that

${\mathbb{E}}[M_{ij}M_{kl}]=\frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})+c\delta_{ij}\delta_{kl},$
(5.1)

where $\delta_{ij}$ is the Kronecker delta function. It is evident that ${\rm
GOI}(c)$ becomes a GOE matrix if $c=0$. The following lemma in [7] provides
the density of the ordered eigenvalues for a GOI matrix.

######  Lemma 5.1.

Let $M$ be an $N\times N$ nondegenerate ${\rm GOI}(c)$ matrix ($c>-1/N$). Then
the density of the ordered eigenvalues $\lambda_{1}\leq\ldots\leq\lambda_{N}$
of $M$ is given by

$\begin{split}f_{c}(\lambda_{1},\ldots,\lambda_{N})&=\frac{1}{K_{N}\sqrt{1+Nc}}\exp\left\\{-\frac{1}{2}\sum_{i=1}^{N}\lambda_{i}^{2}+\frac{c}{2(1+Nc)}\left(\sum_{i=1}^{N}\lambda_{i}\right)^{2}\right\\}\\\
&\quad\times\prod_{1\leq i<j\leq
N}|\lambda_{i}-\lambda_{j}|\mathbbm{1}_{\\{\lambda_{1}\leq\ldots\leq\lambda_{N}\\}},\end{split}$
(5.2)

where $K_{N}=2^{N/2}\prod_{i=1}^{N}\Gamma\left(\frac{i}{2}\right)$.

We use the notation ${\mathbb{E}}_{{\rm GOI}(c)}^{N}$ to represent the
expectation under the ${\rm GOI}(c)$ density (5.2), i.e., for a measurable
function $g$,

${\mathbb{E}}_{{\rm
GOI}(c)}^{N}[g(\lambda_{1},\ldots,\lambda_{N})]=\int_{{\mathbb{R}}^{N}}g(\lambda_{1},\ldots,\lambda_{N})f_{c}(\lambda_{1},\ldots,\lambda_{N})d\lambda_{1}\cdots
d\lambda_{N}.$ (5.3)

###  5.2 Peak Height Distribution of Anisotropic Gaussian Random Fields

Let $A$ be a nondegenerate $N\times N$ matrix, that is ${\rm det}(A)\neq 0$. A
Gaussian random field $\\{X(t),t\in{\mathbb{R}}^{N}\\}$ is called anisotropic
if there exists a function $\varphi:[0,\infty)\to{\mathbb{R}}$ such that the
covariance of $X$ has the form

${\mathbb{E}}[X(t)X(s)]=\varphi\left(\|A(t-s)\|^{2}\right),\quad
t,s\in{\mathbb{R}}^{N}.$ (5.4)

Note that $X$ is stationary but not isotropic if $A\neq I_{N}$. A typical
example is that
${\mathbb{E}}[X(t)X(s)]=\exp\\{-\sum_{i=1}^{N}a_{i}^{2}(t_{i}-s_{i})^{2}\\}$
with $A={\rm diag}(a_{1},\ldots,a_{N})$ and
$\varphi(\theta)=\exp\\{-\theta\\}$, where $a_{i}\neq a_{j}$ for $i\neq j$ and
$\theta\geq 0$. Following the notation in [6, 7], let

$\kappa=-\varphi^{\prime}(0)/\sqrt{\varphi^{\prime\prime}(0)}.$ (5.5)

######  Theorem 5.2.

Let $\\{X(t),t\in{\mathbb{R}}^{N}\\}$ be an anisotropic Gaussian field
satisfying (5.4). Then the peak height distribution of $X$ is given by

$\begin{split}F(u)=\begin{cases}\frac{\int_{u}^{\infty}\phi(x){\mathbb{E}}_{{\rm
GOI}((1-\kappa^{2})/2)}^{N}\left[\prod_{j=1}^{N}\big{|}\lambda_{j}-\kappa
x/\sqrt{2}\big{|}\mathbbm{1}_{\\{\lambda_{N}<\kappa
x/\sqrt{2}\\}}\right]dx}{{\mathbb{E}}_{{\rm
GOI}(1/2)}^{N}\left[\prod_{j=1}^{N}|\lambda_{j}|\mathbbm{1}_{\\{\lambda_{N}<0\\}}\right]},&\quad\text{\rm
if }\kappa^{2}<(N+2)/N,\\\ \frac{{\mathbb{E}}_{{\rm
GOI}(1/2)}^{N}\left[\prod_{j=1}^{N}|\lambda_{j}|\mathbbm{1}_{\\{\lambda_{N}<0\\}}\mathbbm{1}_{\left\\{\sum_{j=1}^{N}\lambda_{j}/N\leq-\sqrt{(N+2)/(2N)}u\right\\}}\right]}{{\mathbb{E}}_{{\rm
GOI}(1/2)}^{N}\left[\prod_{j=1}^{N}|\lambda_{j}|\mathbbm{1}_{\\{\lambda_{N}<0\\}}\right]},&\quad\text{\rm
if }\kappa^{2}=(N+2)/N,\end{cases}\end{split}$ (5.6)

where the expectation ${\mathbb{E}}_{{\rm GOI}(c)}^{N}$ is defined in (5.3).

###### Proof.

Let $Z$ be an isotropic Gaussian field with covariance

${\mathbb{E}}[Z(t)Z(s)]=\varphi\left(\|t-s\|^{2}\right).$ (5.7)

Then

${\mathbb{E}}[Z(At)Z(As)]=\varphi\left(\|A(t-s)\|^{2}\right)={\mathbb{E}}[X(t)X(s)],\quad
t,s\in{\mathbb{R}}^{N},$

implying that $\\{X(t),t\in{\mathbb{R}}^{N}\\}$ and
$\\{Z(At),t\in{\mathbb{R}}^{N}\\}$ have the same covariance structure and
hence the same peak height distribution. We can write $X(t)=Z(At)$, and hence
$\nabla^{2}X(t)=A^{T}\nabla^{2}Z(At)A$. Note that ${\rm
det}(A^{T}\nabla^{2}Z(At)A)=({\rm det}(A))^{2}{\rm det}(\nabla^{2}Z(At))$ and
that $A^{T}\nabla^{2}Z(At)A\prec 0$ if and only if $\nabla^{2}Z(At)\prec 0$.
Therefore, the peak height distribution of $X$ is given by

$\begin{split}F(u)&=\frac{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{X(t)>u\\}}\mathbbm{1}_{\\{\nabla^{2}X(t)\prec
0\\}}\right]}{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}X(t))|\mathbbm{1}_{\\{\nabla^{2}X(t)\prec 0\\}}\right]}\\\
&=\frac{{\mathbb{E}}\left[|{\rm
det}(A^{T}\nabla^{2}Z(At)A)|\mathbbm{1}_{\\{Z(At)>u\\}}\mathbbm{1}_{\\{A^{T}\nabla^{2}Z(At)A\prec
0\\}}\right]}{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}A^{T}Z(At)A)|\mathbbm{1}_{\\{A^{T}\nabla^{2}Z(At)A\prec
0\\}}\right]}\\\ &=\frac{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}Z(At))|\mathbbm{1}_{\\{Z(At)>u\\}}\mathbbm{1}_{\\{\nabla^{2}Z(At)\prec
0\\}}\right]}{{\mathbb{E}}\left[|{\rm
det}(\nabla^{2}Z(At))|\mathbbm{1}_{\\{\nabla^{2}Z(At)\prec
0\\}}\right]},\end{split}$

which is exactly the peak height distribution of the isotropic Gaussian field
$Z$. The desired peak height distribution (5.6) then follows from Corollaries
3.6 and 3.14 in [7] as well as the fact $\kappa^{2}\leq(N+2)/N$ therein. ∎

### Acknowledgments

The author was supported by NSF grant DMS-1902432 and Simons Foundation
Collaboration Grant 854127. He is grateful to Professor Armin Schwartzman from
UC San Diego for helpful discussions and suggestions.

## References

  * [1] Adler, R. J. and Taylor, J. E. (2007). Random Fields and Geometry. Springer, New York.
  * [2] Allard, D., Senoussi, R. and Porcu, E. (2016) Anisotropy Models for Spatial Data. Mathematical Geosciences. 48, 305–328.
  * [3] Bardeen, J. M., Bond, J. R., Kaiser, N. and Szalay A. S. (1985). The statistics of peaks of Gaussian random fields. Astrophys. J., 304, 15–61.
  * [4] Cheng, D., Cammarota, V., Fantaye, Y., Marinucci, D. and Schwartzman, A. (2020). Multiple testing of local maxima for detection of peaks on the (celestial) sphere. Bernoulli, 26, 31–60.
  * [5] Cheng, D., He, Z. and Schwartzman, A. (2020). Multiple testing of local extrema for detection of change points. Electron. J. Statist., 14, 3705–3729.
  * [6] Cheng, D. and Schwartzman, A. (2015). Distribution of the height of local maxima of Gaussian random fields. Extremes, 18, 213–240.
  * [7] Cheng, D. and Schwartzman, A. (2018). Expected number and height distribution of critical points of smooth isotropic Gaussian random fields. Bernoulli, 24, 3422–3446.
  * [8] Cheng, D. and Schwartzman, A. (2017). Multiple testing of local maxima for detection of peaks in random fields. Ann. Stat., 45, 529–556.
  * [9] Cramér, H. and Leadbetter, M. R. (1967). Stationary and Related Stochastic Processes: Sample Function Properties and Their Applications. Wiley, New York.
  * [10] Mallows, C. L. (1961). Latent vectors of random symmetric matrices. Biometrika, 48, 133–149.
  * [11] Pantazis, D., Nichols, T. E., Baillet, S., and Leahy, R. M. (2005). A comparison of random field theory and permutation methods for the statistical analysis of MEG data. Neuroimage, 25(2), 383–394.
  * [12] Taylor, J. E. and Worsley, K. J. (2007). Detecting sparse signals in random fields, with an application to brain mapping. _J. Am. Statist. Assoc._ , 102, 913–928.
  * [13] Van Hecke, W., Leemans, A., De Backer, S., Jeurissen, B., Parizel, P. M., and Sijbers, J. (2010). Comparing isotropic and anisotropic smoothing for voxel‐based DTI analyses: A simulation study. Human brain mapping, 31(1), 98-114.

> Dan Cheng
> School of Mathematical and Statistical Sciences
> Arizona State University
> 900 S. Palm Walk
> Tempe, AZ 85281, U.S.A.
> E-mail<EMAIL_ADDRESS>