Waku: A Family of Modular P2P Protocols For Secure&Censorship-Resistant Communication

Waku is a family of modular protocols that enable secure, censorship-resistant, and anonymous peer-to-peer communication. Waku protocols provide capabilities that make them suitable to run in resource-restricted environments e.g., mobile devices and web browsers. Such capabilities include (i) retrieving historical messaging for mostly-offline devices (ii) adaptive nodes; allowing for heterogeneous nodes to contribute to the network (iii) preserving bandwidth usage for resource-restricted devices, (iv) minimizing connectivity requirements for devices with a limited connection, and (v) enabling efficient, private, economic spam protection for heterogeneous nodes. Waku's modular design and resource-efficient protocols make it superior to its predecessor i.e., Whisper. In this paper, we give an overview of the Waku protocols stack, its architecture, and protocols interaction along with a sample demo scenario on configuring and running a Waku node using nwaku i.e., Waku client written in Nim.

I. INTRODUCTION
Waku is a family of modular protocols designed to enable privacy-preserving peer-to-peer communication for decentralized applications. It is the successor of Whisper [1], the messaging layer of the Ethereum blockchain P2P protocol suite, and outperforms it in scalability and resource efficiency. Below we highlight Waku's features 1 . At the time of this article, Waku is deployed by Status and WalletConnect v2. Peer-to-peer & Censorship-resistant: Waku protocols follow a p2p design where an intended functionality is realized by the resources supplied by independently run peers. Being p2p makes Waku protocols resistant to censorship and a single point of failure also allows utilization of shared infrastructure. Privacy-Preserving: Unlike centralized messaging systems, Waku does not rely on peers' Personally Identifiable Information (PII) e.g., phone number or email address. Also, protocol messages are metadata/PII protected to enable user anonymity and various forms of unlinkability [2]. Useful for generalized messaging: Waku is designed to work in various messaging scenarios including human-to-human e.g., a messenger, or machine-to-machine e.g., state channels. Modularity: Waku protocols are modular hence nodes can adaptively choose protocols to support. We call this concept adaptive nodes explained in 30/ADAPTIVE-NODES RFC [3]. This allows applications to make trade-offs based on the properties they and their users value, e.g., resource usage vs metadata protection; providing useful services to the network vs being just a consumer. Runs anywhere: Waku protocols can be run in environments with limited resources e.g., bandwidth, CPU, memory, disk, and battery, and connectivity (i.e., not being publicly connectable, only being intermittently connected, and mostlyoffline).
II. SOFTWARE ARCHITECTURE We present four core protocols of the Waku stack as illustrated in Figure 1. We note that this is not an exhaustive list. Protocols are referred to by their RFC IDs as in their reference specifications [3] and are mostly implemented over libp2p and associated with libp2p protocol identifiers. Protocols fall into one of the three separate network interaction domains i.e., (a) Gossip domain, (b) Discovery domain, and (c) Request/Reply (Req/Rep) domain as explained next.
Waku Message: Waku protocols are used to communicate a data unit called a Waku Message i.e., 14/WAKU2-MESSAGE [3]. Each message is associated with a content topic used for content filtering. The content topic has no impact on the message routing, but enables fine-grain queries in Req/Reply protocols and can be used to specify an application domain. Gossipsub Domain: Waku uses a publisher-subscriber messaging model to disseminate Waku messages. This is done through the 11/WAKU2-RELAY protocol which is a thin layer on top of libp2p GossipSub [4]. It is a Gossip-based routing algorithm where messages are published to topics and delivered by gossiping to peers subscribed to that pubsub topic. Waku also features an experimental privacy-preserving economic and efficient spam protection mechanism on top of 11/WAKU2-RELAY called 17/WAKU2-RLNRELAY [3]. Discovery Domain Waku supports multiple discovery methods by which nodes become aware of each other's existence. That includes DNS-based discovery as per EIP-1459 [5] by which a new node can retrieve an authenticated, updateable list of peers to bootstrap connection to the network. This method merits high scalability and resource efficiency, though, is centralized hence prone to attacks e.g., censorship. It may be used in conjunction with other ambient peer discovery methods e.g., Node Discovery v5 [6] which is a DHT-based mechanism and, unlike the DNS-based approach, is decentralized. Nevertheless, its high resource usage e.g., CPU and memory must be carefully considered when used for resource-limited devices. Waku also supports the basic form of discovery where a set of static nodes may be specified. Request/Reply Domain: Waku supplies a set of Req/Rep protocols primarily used to get Waku to run in resource-restricted environments e.g., with limited bandwidth or connectivity. Req/Rep protocols have two modes 1) full mode; peer can act as both a requester and a responder 2) light mode; peer is only a requester. The prerequisite of running in full mode is to have support for 11/WAKU2-RELAY protocol.
• 13/WAKU2-STORE provides support for historical messages where it enables querying of Waku messages received through 11/WAKU2-RELAY protocol and stored by other nodes. This is especially useful for mostly offline devices. Nodes running this protocol in full mode are expected to have high availability and ample storage. Waku also provides an experimental fault-tolerant addition to the store protocol called 21/WAKU2-FT-STORE [3] that relaxes the high availability requirement. • 12/WAKU2-FILTER is a lightweight version of 11/WAKU2-RELAY designed for bandwidth restricted devices in which nodes running in light mode can subscribe to nodes running in full mode to only receive Waku messages they desire e.g., based on the pubsub topic or content topic. • 19/WAKU2-LIGHTPUSH is used for nodes with short connection windows and limited bandwidth to publish messages. A receiver of light push message publishes the encapsulated Waku message on the requested pubsub topic. A sample interaction of Waku protocols is visualized in Figure 2. There are three nodes denoted by A-C. The PubSub topic pub1 is used for routing and B is subscribed to that topic. B also supports the 13/WAKU2-STORE protocol. A creates msg1 with content topic content1 and publishes it on pub1. As B is subscribed to that topic, it picks up the message and saves it for possible later retrieval by other nodes. At a later time, C comes online and requests messages matching pub1 and content1 from B. B responds with messages meeting this, and possibly other, criteria.
III. SAMPLE DEMO SCENARIO nwaku [7] is an open-source and reference implementation of the Waku client written in Nim programming language. There are also other implementations of Waku available that support a subset of Waku protocols, namely, js-waku [8] and go-waku [9], written in Javascript and Go, respectively.
Listing 1 illustrates a sample configuration of a Waku node using nwaku client via a set of command-line flags. wakunode2 is the binary executable of nwaku which can be built as instructed in [7]. The node is configured to statically connect to a relay node whose libp2p multi-address [10] is passed in staticnode option. The relay flag enables 11/WAKU2-RELAY protocol over the pubsub topic supplied by topics flag which is /waku/2/default-waku/proto. Similarly, the 13/WAKU2-STORE protocol is enabled by passing the store option. The boolean flag persist-messages is used to distinguish between the full-mode and light mode of the store protocol (true indicates full mode). The store-capacity specifies a limit on the number of persisted Waku messages. The filter and lightpush flags are used to mount the 12/WAKU2-FILTER protocol and 19/WAKU2-LIGHTPUSH in full mode. A Waku node may optionally expose a set of JSON RPC APIs using the rpc flag depending on its supported protocols. More details on these APIs are given in 16/WAKU2-RPC RFC [3]. The corresponding IP address and port can be set using rpcAddress and rpc-port options.


A unified framework on defining depth for point process using function smoothing

The notion of statistical depth has been extensively studied in multivariate and functional data over the past few decades. In contrast, the depth on temporal point process is still under-explored. The problem is challenging because a point process has two types of randomness: 1) the number of events in a process, and 2) the distribution of these events. Recent studies proposed depths in a weighted product of two terms, describing the above two types of randomness, respectively. In this paper, we propose to unify these two randomnesses under one framework by a smoothing procedure. Basically, we transform the point process observations into functions using conventional kernel smoothing methods, and then adopt the well-known functional $h$-depth and its modified, center-based, version to describe the center-outward rank in the original data. To do so, we define a proper metric on the point processes with smoothed functions. We then propose an efficient algorithm to estimated the defined"center". We further explore the mathematical properties of the newly defined depths and study asymptotics. Simulation results show that the proposed depths can properly rank the point process observations. Finally, we demonstrate the new method in a classification task using a real neuronal spike train dataset.

Introduction
The statistical depth is a method to indicate the centrality of a data point with respect to a data cloud. It can provide a center-outward structure that can be used to understand the empirical distribution, similar to the empirical quantiles in univariate data. The concept of depth was firstly proposed by (Tukey 1975) to handle multivariate data. Since then the notion of depth has been widely studied by mathematicians and statisticians. Different forms of depth were proposed and analyzed, based on different usages and criteria. The types of data that those depths were mainly applied to were multivariate data and functional data.
Some well-known multivariate depths include the half-space depth (Tukey 1975), the convex hull peeling depth (Barnett 1976), the Oja depth (Oja 1983), the simplicial depth (Liu et al. 1990), and the Mahalanobis depth (Liu & Singh 1993). A study by Zuo & Serfling (2000) discussed desirable properties of multivariate depth, which include affine invariance, maximality at the center, monotonicity relative to the deepest points, and vanishing at infinity.
As a generalization of multivariate data from finite dimension to infinite dimension, functional data received a lot of attentions recently and functional depths were also well studied.
To name a few, Cuevas et al. (2007) proposed the h-depth, which applied the Gaussian kernel and L 2 norm to form the depth. Cuesta-Albertos & Nieto-Reyes (2008) extended the idea of half-space depth and defined the random Tukey depth. In 2009, the band depth and modified band depth were introduced by López-Pintado & Romo (2009), which were very commonly used in functional data problems. Similar to the work by Zuo & Serfling (2000), Nieto-Reyes et al. (2016) examined the desirable properties of functional depth, including distance invariance, maximality at center, strictly decreasing with respect to the deepest point, upper semi-continuity, and receptivity to convex hull width across the domain and continuity.
As a special type of data, observation from an orderly temporal point process is made up by an ascending sequence of event times. Such observation contains two types of randomness: the number of events and the time locations of these events. If we treat each observation as a vector of the event times, then the dimension of this vector will be a random variable. Given this dimension, this vector will also be random with ascending entries.
We point out that the study on the notion of depth of point process is relatively new. So far the only existing work in this area was done in Liu et al. (2017), where the depth structure was built in two steps: For an observation s, 1) estimating the probability of getting the number of events P (|s|); 2) given the number of events, compute the conditional depth D(s | |s|). The estimation of probability was done through a normalized Poisson mass function and the conditional depth adopted the Mahalanobis depth for multivariate data. The final depth was the multiplication of these two with a weight power r: D(s) = P (|s|) r D(s | |s|). This depth structure involves both types of randomness and satisfies good mathematical properties such as invariant to time-shift and linear transformation, monotone on rays, and upper semi-continuous. However, the combination of the randomness on the number of events and the event time distribution is not natural that the impact of the number of events on the final depth should be adjusted by the hyper-parameter r. The selection of the hyper-parameter could be tricky because an inappropriate value will make one of the two randomnesses dominates the depth value. In addition, the two-step procedure deals with the number of events and the event time independently. A more desirable method should be able to combine both steps in one framework, where the two types of randomness can be measured at the same time.
To deal with these issues, we propose another approach to process the point process observation: a "transformation" through smoothing method (Wand & Jones 1994). The idea is to smooth the observed point process sequence with an appropriate kernel function.
Through this operation, the event time vector will be replaced by a function curve. In Section 2.1.2, we will see that the point process observation and the function curve are oneto-one matched if the kernel function used for smoothing satisfies certain mild conditions. According to this bijective relation, we are able to apply the methods on functional data to the point process observations, such as metrics on functions and functional depths. In this way, both types of randomness are taken into account under one framework. For individual event time, it determines the location of the kernel function curve. Because we take the sum of the kernel functions in the smoothing procedure, a larger number of events will enhance the smoothed curve vertically. In addition, irrespective of the number of entries in the vector, its smoothed curve will be just one function. That is, vectors with different dimensions can be natuarally compared. Therefore, this approach will not suffer the same issues in the framework proposed by Liu et al. (2017) and it has the advantage to utilize existing depth methods for functional data.
The rest of this manuscript is organized as follows. In section 2, we will introduce the details of transforming a vector to a smoothed curve. We will also propose a proper metric for the space of point process observations, followed by a discussion of the properties for this metric. With this metric, we will continue to define new depth methods for the point process observations. In section 3, we will examine the asymptotic theory of the proposed depths. In section 4, we will apply our methods on simulations and experimental datasets to validate the effectiveness of the new method. Finally in section 5, a summary of the paper will be shown, followed by the future work. All mathematical proofs and algorithmic details are given in appendices in the supplementary material.

Equivalent Representation via Function Smoothing
Our goal is to rank point process observations on a finite interval [0, T ] via a kernel smoothing method. We will at first introduce the type of smoothing kernels for this purpose.

Kernel Functions
The basic idea of kernel smoothing on an orderly temporal point process is to assign a kernel to each observed point and then sum over all the assigned kernels to get a smooth function. Often the kernel is a probability density function of a given distribution such as a Gaussian kernel. The kernel function will depend on the time interval [0, T ] for the point process and we denote it as K(·; T ). In general, we propose to use any kernel function which satisfies the following four basic requirements. These requirements are needed in order to achieve good mathematical properties for the new depth function: 1. Continuous and non-negative: 2. Positive at zero: K(0; T ) > 0; 3. Linear independence with shifting: for any n (n ∈ N) different values t 1 < t 2 < · · · < t n , we have: = 0 for any x ∈ [0, T ] ⇐⇒ α 1 = · · · = α n = 0; 4. Scale invariance: for any x ∈ [0, T ] and constant α > 0, K(αx; αT ) ≡ K(x; T ); We define a kernel function to be "proper" if it satisfies the above four conditions. To have an example of a proper function, we can consider the following Gaussian kernel function: where c 1 , c 2 are two positive constants. This result is given in the following lemma, where the proof is given in Appendix A.
Lemma 1. The Gaussian kernel in Eqn. (1) is a proper kernel.

Proper Metric
Before moving on to the depth function, we look for a proper metric to measure the difference between two observed point processes. Since the two point processes, in general, may have different number of events, multivariate measurements cannot be applied directly due to inconsistency in dimensions. As smoothing functions are included in this depth, here we will use the distance on smoothed processes as the distance on point processes. Let the time interval of the point process be [0, T ]. We define the set of observed point process with given dimension l > 0 as: l=0 Ω l is the space of all point processes. Each process x = (x 1 , x 2 , · · · , x l ) ∈ Ω l with l > 0 can be represented using a Dirac delta function in the form: x = l i=1 δ(· − x i ). Let K(·; T ) be the smoothing kernel. Then the smoothed process is a function on [0, T ] in the form: In general, the space of smoothed processes for l events with l > 0 is We point out that the correspondence between a point process and its smoothed version is one-to-one. This is given in the following lemma (the proof is in Appendix B).
Lemma 2. The smoothing procedure is a bijective mapping from Ω to F. Now we are ready to define a metric on the point process space. This definition is based on the classical L p norm on the smoothed processes. This is formally given as follows: Definition 1. For any two point processes s, t ∈ Ω on [0, T ] and the correspondent f s , f t ∈ F given by Eqn.
(2), we define a distance function d K,p between s and t as: Note that the smoothing kernel K(·; T ) and L p norm can influence the distance value.
However, we can prove that the distance d K,p in Eqn.
(3) is a proper metric on Ω for any K and p. This is given in the following theorem (the proof is in Appendix C).
Theorem 1. If the smoothing kernel K satisfies the four conditions to be proper, then the function d K,p in Definition 1 is a proper metric in the point process space Ω. That is, it satisfies non-negativity, identity of indiscernibles, symmetry, and triangle inequality.
Based on Lemma 2 and Theorem 1, we know that the point process space Ω and smoothed process space F are one-to-one, and an L p norm on F can be used to define a proper metric on Ω. This metric helps solve the problem in Ω where a conventional vector metric cannot be directly used as different processes may have different cardinalities. Based on this result, we need to answer two natural questions: 1. For two point processes with the same cardinality, if the event times in the two processes are very close, will the distance d K,p between them be close, too?
2. Conversely, if the distance d K,p between two point processes is close, will they have the same cardinality and their events times are also close?
Question 1 examines if the d K,p distance is continuous with respect to the event times.
We claim that this is true and the conclusion is stated in Proposition 1 as follows. The detailed proof is given in Appendix D.
Proposition 1. For any k ∈ N, suppose y = (y 1 , y 2 , · · · , y k ) is an observed point process k ), n = 1, 2, · · · be a sequence of processes in Ω k . If lim n→∞ x (n) = y, or equivalently, lim n→∞ x (n) i = y i , i = 1, 2, · · · , k, then for any observed point process z = (z 1 , z 2 , · · · , z l ) ∈ Ω, Question 2 examines the inverse continuity. We claim that this inverse continuity is also true and the result is stated in Proposition 2 (see Appendix E for a detailed proof).

Desirable Properties for Depth on Point Process
In general, a depth definition is defined to provide a measure of centrality of a given data point within a data cloud. To measure if such goal is achieved, one often examines desirable mathematical properties corresponding to the centrality measurement. For instance, Zuo & Serfling (2000) proposed the desirable properties for multivariate depth. Later, Nieto-Reyes et al. (2016) provided the desirable properties for functional depth. The commonly studied properties are (1) linear invariance, (2) vanishing at infinity, (3) maximality at the center, and (4) monotonicity.
Motivated by previous studies on desirable properties for multivariate and functional depths, we propose the following five desirable properties for the depth on point process: Suppose D(s; P S ) is the depth function for observed event time vector s with respect to the probability space (Ω, F , P S ) for a random point process S on interval [0, T ], we expect that D satisfies P1 Continuity: D(s; P S ) is continuous with respect to s in Ω.

P2 Linear invariance: For any event time vector
P3 Vanishing at infinity: When the number of events goes to infinity, the depth should go to 0: D(s; P S ) → 0 as |s| → ∞.
P4 Unique maximum at the center: There exists s c ∈ Ω, such that D(s c ; P S ) = max s∈Ω D(s, P S ).
Also, for any x( = s c ) ∈ Ω, D(x; P S ) < D(s c ; P S ). We refer to s c as the "center".
P5 Monotone decreasing from the center: For any s 1 .s 2 ∈ Ω, suppose the center is s c , if We point out that the "center" in P4 was originally defined as a point with symmetry in the data cloud. For example, when dealing with zero-mean multivariate normal samples, we will consider the origin as center, and the corresponding Mahalanobis depth is uniquely maximized at this point. However, for point process data in a finite domain [0, T ], it is difficult to define a geometrically symmetric central point in Ω. Notwithstanding, we expect the proposed depth has a unique maximum point, which results in our notion of "center".
Here we examine if the above five properties are satisfied by the generalized Mahalanobis depth on point process (Liu et al. 2017): 1) If two point processes are close to each other, they should have the same number of events and close event times in the time order. As the classical Mahalanobis depth is continuous with respect to the input vector, the generalized Mahalanobis depth on point process satisfies P1; 2) Linear transformation on the time will not change the number of events, so P (|s|) r will keep the same. Because the classical Mahalanobis depth is linear invariant, P2 is also satisfied; 3) The Poisson mass function will go to 0 as the input goes to infinity, so P (|s|) → 0 as |s| → ∞. Thus, P3 holds, too; 4) D(s | |s|) has maximum value to be 1. However, this maximum in general can be achieved at multiple processes. For example, if the total intensity Λ is an integer, both the population mean with dimension Λ and the population mean with dimension Λ − 1 will achieve the maximum depth. Thus, P4 does not hold as the solution to the maximum is not unique; 5) As a result, P5 does not hold either. In summary, the generalized Mahalanobis depth on point process only satisfies P1 to P3.
Similar to the discussion about the Mahalanobis depth, we will refer to P1 -P5 to explore these properties for other depth functions in the following sections of this paper.

h-depth on Point Process
We have defined a proper metric in the point process space Ω. This metric is based on a smoothing procedure where a point process is equivalently represented by a function in F.
In this section, we will exploit a commonly used functional depth, called h-depth, to define an "h-depth" for point processes in Ω.
At first, we review the notion of h-depth (Cuevas et al. 2007) for functional data. Assume X is a functional random variable in the probability space (Λ, F Λ , P Λ ), where Λ is a subset of L 2 ([0, T ]). Let · denote the conventional L 2 norm. For any function z ∈ Λ, its h-depth is defined as: where G h (·) = exp(− t 2 2h ) is a modified Gaussian kernel with parameter h (> 0) and X is a random function on the probability space (Λ, F Λ , P Λ ). To simplify notation, we use HD(z) to represent HD(z; P Λ ) when the measure P Λ is implicitly known. Given a set of i.i.d.
random functions X 1 , X 2 , · · · , X n ∈ Λ, the sample version of the h-depth of z is given as ).
One important property of the h-depth in Eqn. (4) is its continuity. This is given in the following lemma, where the proof can be found in Appendix F.
Lemma 3. If the functional random variable X satisfies E X < ∞, then HD(·) is a continuous function on Λ.
Based on the notion of h-depth on functional data (in Eqn. (4)) and the proper metric on point process (in Eqn. (3)), we can formally introduce the h-depth to point process observations as follows: Definition 2. Let S be a random point process on [0, T ] in the probability space (Ω, F , P S ).
The h-depth of any s ∈ Ω is defined as: where f s and f S are smoothed curves for s and S (by Eqn. (3)) and G h (t) = exp(− t 2 2h ).
Definition 2 provides the population version of the depth. In practice, we should follow the sample version of h-depth on the given observations, which is given below.
be a sample of event time vectors from a point process on [0, T ]. The empirical h-depth of any s ∈ Ω is defined as: where f s and f S i are smoothed curves for s and With the definitions for both the population and empirical versions of the h-depth on point process, we now examine the mathematical properties based on the discussion in section 2.2. The basic result can be summarized in the following proposition. • P3: D(s; P S ) → 0 when |s| → ∞.
The detailed proof is given in Appendix G. Note that P4 and P5 are not satisfied because h-depth in general may not have a unique maximum point in the data cloud.

Modified h-depth on Point Process
The center has been a critical notion in statistical depths. However, as we have pointed out in Sec. 2.3, the h-depth for a point process may not have a center. In this section, we propose to modify the h-depth by including a "center"-based process in the definition.

Definition and Properties
Many commonly used depth functions have the "center" as the point with maximum depth.
Based on this idea, we propose a center-based new depth function on point process using the smoothing method. The formal definition is given as follows: Definition 4. Let s c be a given "center" point process in Ω on [0, T ]. For any s ∈ Ω, its center-based h-depth is defined to be: where h > 0 is a parameter. f s and f sc are smoothed processes of s and s c , respectively. )], with X to be a functional random variable. If the order of the exponential function and the expectation is switched, then we have ). In this way we obtain the center-based h-depth D(s; s c ).
Remark 2: The L 2 norm in the definition can be generalized to L p , with 1 ≤ p < ∞.
With a center process given in the definition, this modified h-depth for point process is expected to satisfy more desirable mathematical properties than the classical one. Indeed, all the desirable properties in section 2.2 are satisfied and the details are given in the following proposition. • P5: For any s 1 , s 2 ∈ Ω, if d K,2 (s 1 , s c ) < d K,2 (s 2 , s c ), then D(s 1 ; s c ) > D(s 2 ; s c ). The detailed proof for each property is provided in Appendix H. Note that this is a clear advantage over the classical h-depth method, where only properties P1-P3 are satisfied.

Estimation of the Center
In practice, there is no prior knowledge of the center of point process, so a proper estimation of the center will be needed. In this study, we adopt the common notion of the Karcher mean in a metric space (Grove & Karcher 1973). Since we have a well-defined metric between any two point processes in Definition 1, the Karcher mean can be defined as: Definition 5. Let S be a random point process on [0, T ] in the probability space (Ω, F , P S ).
The Karcher mean of S is defined to be where d K,2 (·, ·) is the metric with kernel function K(·; T ) and L 2 norm in Definition 1.
Note that, in general the solution of the minimization may not be unique, so µ K is a set of points. Based on the Karcher mean in Definition 5, its empirical version is given as follows: The empirical Karcher mean also may not be unique, and thereforeS K is a set of solution points, too. For simplicity, we define the following notation: where SSD stands for "Sum of Squared Distance" and t ∈ Ω is the event time vector. With this definition, the empirical Karcher mean is justS (N ) K = arg min t∈Ω SSD(t; S 1 , S 2 , · · · , S N ). For any element inS We will introduce a theoretically-proven convergent algorithm to estimateS K provides an effective searching range. Such algorithm is a combination of simulation annealing and line search. We will at first briefly review the simulation annealing and line search, and then introduce an approach to combine these two methods.
• RJMCMC annealing Simulated annealing (Geman & Geman 1984, Van Laarhoven & Aarts 1987) is a method that makes optimization feasible if we are able to generate samples. To find the global minimum of function f (x), we can build a density function π where T i is a pre-determined decreasing sequence of positive numbers called temperature and lim i→∞ T i = 0 (e.g. T i = C/ log(1 + i) for constant C > 0). Under weak regularity assumptions on f , the density π i will be concentrated on the global minimum points of f . After a large number of iterations of computing π i and sampling from it, we will have samples passing the global minimum point.
Simulation annealing provides a method to find the minimum of a function, but to apply it on the optimization of SSD, a tool that is able to generate vectors with different dimensions will be needed. In this paper, we adopt the reversible-jump Markov Chain Monte Carlo (RJMCMC) method, which was first introduced by Green (1995). Its main idea is the Metropolis Hasting method (Metropolis et al. 1953), but generalizes to allow the candidate to have different dimensions. In this way, it can be used to get samples from densities where the dimension of the random vector is not fixed. One condition to use the RJMCMC method is that the number of possible dimensions should be finite.
According to Theorem 2, the range of dimension is 1, · · · , D (N ) K , so RJMCMC is feasible for the center estimation.
Combining the two methods, we are able to find the Karcher mean through sampling from T i SSD(t; S 1 , S 2 , · · · , S N )]. Because of RJMCMC, the annealing process can search over different dimensions automatically. An example RJMCMC annealing algorithm is provided in Appendix J, where the input to the algorithm can be adjusted before application.
• Line search According to Theorem 2, because the dimensions to search is finite 1, 2 · · · D (N ) K , optimizing within each dimension and comparing the output across dimensions will be another feasible approach, which corresponds to the "line search" idea. An advantage of this method is that the gradient of the SSD can be computed explicitly with the input dimension fixed, as shown in the following proposition.
Proposition 5. If the dimension of the input event time vector t is given and the distance function is d k,2 , then the gradient of SSD with respect to t will be: where t and S i are event time vectors. x is scaler and |t| is the dimension of vector t.
All operations are done element-wise.
If the modified Gaussian kernel K G (t; T ) = c 1 e − c 2 T 2 t 2 is used to smooth the process, then the gradient can be further simplified to: with the function g(·, ·) to be where Φ is the cumulative distribution function for standard Normal distribution N(0, 1); 1 d = [1, 1, · · · , 1] T ∈ R d . All operations are done element-wise.
The proof and computational details are provided in Appendix L. Due to this property, gradient-based methods, such as stochastic gradient descent, can be applied to efficiently do the optimization within a given dimension. Then a comparison of the output will return the solution to the minimum SSD, which is the empirical Karcher mean. An example line search algorithm based on gradient method is provided in Appendix K.
• Combined method Both RJMCMC annealing method and line search method can be used to estimate the empirical Karcher mean, but they have some disadvantages. For the RJMCMC annealing, although it can search over dimensions automatically, it converges slowly to the optimal solution. For line search, it can converge fast to the solution given the dimension, but it needs to search every dimension in the dimension range. As a result, the time cost for the entire line search algorithm is also very large.
Based on their characteristics, if we combine the two methods: do RJMCMC annealing first as a "pre-train" process to locate a range of optimal dimensions and then use line search to do optimization within the narrowed dimension range, the computation efficiency is expected to be significantly improved. Moreover, the output event time vectors from the RJMCMC annealing can be used as the initial values for the line search.
In this way, the number of iterations in the line search will be reduced and it can also help avoid local optimal points. The combined algorithm following this idea is shown in Algorithm 1. The application of this algorithm on simulations and real data will be shown in section 4.

Asymptotic Theory
In this section, we will show that the empirical depth converges to the population depth as the sample size goes to infinity for h-depth and modified h-depth, defined in section 2.
Algorithm 1 Combined method to find empirical Karcher mean Input: the observed event time vectors S 1 , S 2 , · · · , S N .
For RJMCMC annealing: The maximum number of iterations n max ; Initial value x 0 ; Initial dimension k 0 to be the dimension of x 0 ; the upper bound of dimension D For combination: The number of dimensions to be kept in RJMCMC annealing: d r ∈ N.
(1) Pre-train: Use RJMCMC annealing (such as Algorithm 2) with input S 1 , S 2 , · · · , S N , n max , x 0 , k 0 and D is the corresponding sample version with sample size N. Then, for any input s ∈ Ω, This theorem implies the convergence of sample h-depth to the population h-depth, and can be easily proven using the Strong Law of Large Numbers.
For the modified h-depth, the convergence of the depth value will be based on the convergence of the estimated center. As defined in section 2.4.2, the Karcher mean can be treated as the solution in Ω that minimizes the average of squared distances, where the distance is d k,2 with L 2 norm. In fact, we can generalize the Karcher mean definition by using any d k,p distance for p ≥ 1, given in the following form: Definition 7. Let S be a random point process on [0, T ] in the probability space (Ω, F , P S ).
The generalized Karcher mean of S is defined to be where d K,p (·, ·) is the metric with kernel function K(·; T ) and L p norm where p ≥ 1.
It is apparent that Definitions 5 and 6 are special cases with p = 2 in the generalized definitions, respectively. Then the convergence theorem is given as follows: 1. The minimum average of squared distance in {S i } N i=1 will converge to the minimum expected squared distance in (Ω, F , P S ) almost surely, in other words, 2. The empirical generalized Karcher mean converges to the generalized Karcher mean almost surely, in other words,S (N ) The almost-surely convergence from setS is the Kuratowski upper limit (Kuratowski 2014):

Application Results
In this section, we will apply the proposed depth framework, h-depth and modified h-depth, on point process observations in simulations and real experimental data.

Homogeneous Poisson Process
We will at first illustrate the depth methods on observations from HPP(λ) in a finite time interval, where λ is a constant mean value. 100 independent realizations from HPP(0.045) on [0, 100] are generated, where the raster plot of these processes is shown in Figure 1(a).
Basically, the number of events in each process follows a Poisson distribution with mean 4.5 (varying from 1 to 10 in these 100 realizations), and the event time in every observation is uniformly distributed on [0, 100]. We then smooth the realizations using a kernel in Eqn.
Based on the h-depth and modified h-depth in section 2, we compute the sample depth value for each given process where the constant h is set to T = 100. For the modified h-depth, the center is determined in four ways: set by prior information, estimated by RJMCMC annealing, line search and combined method. As the Poisson process is homogeneous, we can intuitively set the center as (20,40,60,80). The output of the center estimation is shown in Table 1. We can see that all the estimated centers have lower SSD than the intuitive center. In particular, the combined method has the superior performance. It has the smallest SSD and also the smallest time cost. We will only use this method for the modified h-depth in comparison with other depth methods.  We then display the corresponding smoothed processes in Figure 3. Note that the generalized Mahalanobis depth is not included as it is not based on smoothing functions. All three plots demonstrate a center-outward decreasing depth value structure that observations whose smoothing curves are closer to the middle of all the curves will have larger depth values. Overall, these simulation results are reasonable and consistent with the basic notion of center-outward ranks.

Inhomogeneous Poisson Process
We will further apply the depth on point process using observations from IPP[λ(t)] to evaluate the perforamnce. The rate function λ(t) is chosen as λ(t) = 3φ(t; 25, 10) + 2φ(t; 75, 10) and the event time interval is still [0, 100], where φ(·; µ, σ) is the density for normal distribution with mean µ and standard deviation σ. As λ(t) contains two peaks, the simulated event times will mainly distribute around 25 and 75. The kernel function is also the modi-  Table 2. It can be seen that the three estimated centers are very similar. The line search method and the combined method both have the lowest SSD, whereas the latter one is much more efficient. That is, the combined method results in the superior performance, and therefore we will only use this method for the modified h-depth in comparison with other depth methods.   (b) By the modified h-depth with center estimated using the combined method.

Method
After computing the depth for all the observations, we obtain the ranking result. Figure   5 shows the top 5 and bottom 5 ranked observations. Similar to the HPP case, we include the result from the generalized Mahalanobis depth with r = 1 for comparison. It can be seen that the outputs from h-depth and modified h-depth are similar -observations with event times around 25 and 75 and number of events around 5 have larger depth. In contrast, for the Mahalanobis depth, we can also observe that processes with number of events around 5 have larger depth. However, rather than concentrating at two peak time locations, the top 5 processes have event time more uniformly distributed. It indicates that h-depth and modified h-depth can better capture the pattern inside the sample. If further taking into account the true λ(t) that the peak at 25 is higher than the peak at 75, we should expect to see more event times concentrated around 25 than around 75. Comparing In summary, based on the result from the HPP and IPP simulations, we can see that the proposed h-depth and modified h-depth both can properly build a center-outward rank on the given point process data.

Experimental data application
We will examine the proposed depth framework using a real spike train recording, where the data were included in the Quantitative Single-Neuron Modeling Competition 2009 (Naud et al. 2009) and were accessible at http://dx.doi.org/10.6080/K0PN93H3. The detailed experiment description can be found in (Carandini et al. 2007).
Briefly, the experiment was performed on rhesus monkeys that retinal input (visual stimulus) was applied and extracellular potentials were recorded for both the retinal (presynaptic) and the geniculate (post-synaptic) simultaneously. There were 10 seconds stimulus, where the first 5 seconds were the same for all trials and the last 5 seconds were unique for each trial. In total, 76 trials were performed. Since the task in the competition was to predict the post-synaptic spikes given the pre-synaptic spikes, only 38 pairs of (post-synaptic, pre-synaptic) were fully given. We will use these spike train data to test the smoothing depth framework through a classification task: In the training set, we have labeled observations for pre-synaptic (abbreviated as "pre-group") and post-synaptic (abbreviated as "post-group"). In the testing set, we compute the depth for each observation in the two groups. The group with a larger depth value will be the predicted label. Before the classification, an exploratory data analysis is performed. The spike times are within the 10 seconds range and the point process plot is shown in Figure 7. 10 example pairs of pre-group, post-group observations are shown. As stated in the experiment description, the stimuli were the same for the first 5 seconds and different for the last 5 seconds. This well explains why the signals have very similar pattern in the first 5 seconds.
Therefore, it is reasonable to separate the trial into two parts, first 5 seconds and last 5 seconds, and we can smooth them in different ways.
After separation, we apply the modified Gaussian kernel to smooth the observations with c 1 = 1, c 2 = 100 for first 5 seconds and c 1 = 1, c 2 = 50 for last 5 seconds. The smoothing curve plot is given in Figure 8. Using the combined method in section 2.4.2, we estimate the Karcher mean for the two groups and the corresponding smoothing curves are shown as thick lines. It can be seen that the Karcher mean curves are able to catch the shape patterns in the first 5 seconds as the observations are nearly aligned, while in the last 5 seconds, the Karcher mean is nearly flat since the observations are very noisy. The classification result is shown in Table 3. All methods are able to provide reasonable classification power on the spike train data. The overall results on the first 5 seconds is better than the last 5 seconds. For the last 5 seconds, it can be seen that the modified h-depth has the best accuracy (76.32%) and F1 scores ( Table 3: The classification result. F1 score (pre) and F1 score (post) are computed according to the recall and precision on pre-group classification and post-group classification.
The boldface indicates best results in each column.
best accuracy (81.58%), while the modified h-depth is better on the post group (F1 score 82.05%) and the other two are better on the pre-group (F1 score 81.58%). In summary, we can see that the smoothing depth framework provides relatively accurate classification in the practical point process observations.

Summary and Future Work
In this paper, we at first proposed a kernel smoothing representation for point process observations. We found that the smoothing procedure builds a bijective mapping between them. We then defined a proper metric distance bewteen the smoothing curves and transformed the problems from point processes to smooth functions. Based on the notion of h-depth on functions, we proposed two methods to form the depth structure on point process observations, i.e. the classical h-depth and the modified h-depth with center given.
We then proposed a center-estimation method for the modified h-depth, with the idea of minimizing the Karcher mean through a combination of the RJMCMC annealing and line search. During the tests on the simulated data, both classical h-depth and modified hdepth resulted in reasonable outcomes: a center-outward decreasing depth structure that observations with smoothing curves closed to the center will have higher depth. In the real neuronal spike train data, we tested the new depth structures with a depth-based classification task and both methods resulted in accurate classifications.
Our investigation is only the starting point of the smoothing curve exploration and there is much to be further studied in the future. We have proposed the modified Gaussian kernel that satisfies the four basic requirements. We will explore other kernel functions that also satisfy the requirements. In the modified h-depth, a combination of RJMCMC annealing and line search has been given as the method to estimate the Karcher mean, which serves as the "center". There should be other ways to define the "center" and other estimation methods. These are all very good topics to further develop depth function on point process. Finally, more functional depth methods, other than h-depth, can be adopted to smoothed point processes, such as the band depth and modified band depth.
Their ranking performance needs thorough investigations.


Benchmarking optimization algorithms for auto-tuning GPU kernels

Recent years have witnessed phenomenal growth in the application, and capabilities of Graphical Processing Units (GPUs) due to their high parallel computation power at relatively low cost. However, writing a computationally efficient GPU program (kernel) is challenging, and generally only certain specific kernel configurations lead to significant increases in performance. Auto-tuning is the process of automatically optimizing software for highly-efficient execution on a target hardware platform. Auto-tuning is particularly useful for GPU programming, as a single kernel requires re-tuning after code changes, for different input data, and for different architectures. However, the discrete, and non-convex nature of the search space creates a challenging optimization problem. In this work, we investigate which algorithm produces the fastest kernels if the time-budget for the tuning task is varied. We conduct a survey by performing experiments on 26 different kernel spaces, from 9 different GPUs, for 16 different evolutionary black-box optimization algorithms. We then analyze these results and introduce a novel metric based on the PageRank centrality concept as a tool for gaining insight into the difficulty of the optimization problem. We demonstrate that our metric correlates strongly with observed tuning performance.

I. INTRODUCTION
G RAPHICS Processing Units (GPUs) have revolutionized the HPC landscape in the past decade [13], and are seen as one of enabling factors in recent breakthroughs in Artificial Intelligence (AI) [23]. GPUs originated as processors for gaming and then adapted to more general workloads as co-processors in many HPC systems. Over the past decade, GPUs have started to again penetrate new markets such as IoT devices [35] and autonomous vehicles [28]. The range of applications of GPUs as such continues to expand. Because of their relatively low cost with respect to their parallel processing power, more and more supercomputers come equipped with GPUs, and in 2020, the majority of modern supercomputers use GPUs [59] as the major source of compute power.
The sections of code that run on a GPU, called kernels, can be challenging to configure such that they run efficiently for a varying combinations of datasets and GPU architectures [63]. The kernel parameters can be split into those defined by the program, and those that are a consequence of the underlying architecture and models behind the GPU. The hardwarespecific parameters define how the thousands of threads in a GPU are grouped. An ineffective layout can cause underutilization of GPU resources. In general, the computational efficiency can drop by an order of magnitude depending on certain implementation choices. Typically, only a small subset of the possible configurations lead to a large increase in performance [53]. Therefore, it is vital to be able to select an efficient kernel configuration.
The search space for this problem is formed by all feasible combinations of GPU kernel parameters. This space is discrete and non-convex [45], making it hard to carry out the optimization. For most GPU kernels used in practice, the size of this search space is such that traversing the options by hand or brute-force is infeasible. An additional complication in optimizing kernel parameters is that evaluating the performance of each configuration requires costly recompilation and test runs. Furthermore, the same GPU kernel often requires re-tuning for different input data, hardware, or after changes to the code [18], [24], [40], [41]. Large throughput pipelines often rely on computationally expensive GPU kernels that consume large amount of resources [54], [55], and cannot be tuned exhaustively due to the aforementioned reasons.
Automatic performance tuning (auto-tuning) techniques rely on empirical results and feedback to optimize the kernel parameters with respect to desired performance metrics. These techniques aim to be widely applicable across architectures. For this reason, auto-tuning can be used to find configurations with increased performance for GPU programs. As the search space for the auto-tuning task depends on various aspects (kernel source, code layout, input data, GPU-architecture), the optimization framework must deal with a broad variety of search spaces and constraints. We, therefore, treat the problem as a black-box optimization task. This raises the question of which optimization algorithm is best suited to find highly efficient settings for GPU kernels, and how these optimization algorithms need to be configured to tune GPU kernels.
The main contribution of this work is to determine which optimization algorithms produce the fastest GPU kernels for different tuning-time ranges. To do so, we conduct a survey of 16 evolutionary optimization algorithms for 9 different NVidia and AMD GPUs, and run 3 real-world applicable benchmark kernels. We select our benchmark problems such that we are able (given ample time) to compute the entire search space, and make these spaces publicly available. To benchmark the GPU kernels we use the Kernel Tuner package [61]. We use arXiv:2210.01465v1 [cs.DC] 4 Oct 2022 the wide range of optimization algorithms present in Kernel Tuner for a large-scale comparison, and provide favourable default hyperparameters for GPU tuning for each algorithm. In addition, we extend Kernel Tuner with several highly-efficient optimization algorithms, including iterative local search (ILS) and dual annealing that cannot be found in any other generic auto-tuning framework.
Secondly, we aim to quantify tuning difficulty for these seemingly challenging and capricious search spaces. To do so, we introduce the fitness flow graph (FFG), which is a network of the points in the search space, with directed edges between neighbours with a better fitness. By computing the likelihood of local search walks terminating in good local minima, we use FFGs to better understand the discrepancies between optimization algorithms, and subsequently tailor them to better suit GPU tuning. In addition, FFGs can help explain the differences across GPU manufactures, architectures, and kernel programs, and such knowledge can help steer future development. To quantify tuning difficulty per kernel, we introduce a novel metric based on Google's PageRank algorithm [3], [43].
This work is structured as follows. In section II we discuss existing GPU kernel tuning approaches. In section III we introduce the preliminaries on GPU kernels, and describe the optimization algorithms that are considered in this survey. In section IV, we describe certain implementation details of Kernel Tuner and our Python optimization package BlooPy. We discuss the setup of our experiments in section V. In section VI we tune the hyperparameters of the algorithms, and present our findings on optimization algorithm performance. In section VII we introduce fitness flow graphs (FFGs) and quantify tuning difficulty for kernel search spaces. Finally, we present our conclusions in section VIII.

A. Automated performance tuning
It is well-known that GPU tuning can yield considerable gains in computational efficiency and utilization for largescale, high-throughput pipelines that run on compute clusters. As an example, we mention the AMBER pipeline [54], [55], which is used to detect Fast Radio Bursts (FRBs) and other single pulse radio transients in astronomy. The pipeline has a throughput of 2 TB/s, and uses a large amount of resources. Benchmarking a single configuration is expensive, and the search space consists of millions of configurations, meaning that sophisticated tuning approaches have to be developed.
Research in automated performance tuning (auto-tuning) can be grouped into two main categories: (1) auto-tuning compiler-generated code optimizations [20], [44], [48], [57], and (2) software auto-tuning [24], [68]. Ashouri et al. [2] wrote an excellent survey on machine-learning methods for compiler-based auto-tuning. In this paper, we limit our scope to (2), i.e., optimizations methods for software auto-tuning, which is sometimes referred to as automated design space exploration [36]. Software auto-tuning allows developers to automatically optimize individual functions and allows, for example, to tune for entirely different implementations and parallelizations that solve the same problem.
A number of generic auto-tuning frameworks have been introduced in recent years. OpenTuner [1] was one of the first generic software auto-tuning frameworks, supporting a number of different search optimization algorithms, but with no support for tuning individual GPU kernels. GPTune [29] and HyperMapper [36] are recently proposed frameworks that both use Bayesian Optimization for auto-tuning on different platforms, but do not target GPUs.
Grauer-Gray et al. [11] have applied auto-tuning to the highlevel directive-based HMPP framework, which can compile to CUDA or OpenCL code. They demonstrate significant performance improvements using auto-tuning over unoptimized HMPP kernels in the PolyBench benchmark suite. Wang et al. [65] take a similar compiler-based approach to automatically convert shared memory OpenMP applications into OpenCL code for GPUs. They use a machine-learning approach, which based on the number of compute operations and memory accesses in the kernels, predicts the best performing hardware platform to execute the kernels either the multi-core CPU using OpenMP, or on the GPU using OpenCL. Hou et al. [16] proposed a data-sensitive auto-tuning framework for sparse matrix vector (SpMV) multiplication that automatically finds the best parallelization strategy. They use a two-step machine learning approach in which they first determine the optimal way to group data into bins and then select the most suitable kernel to process the rows in each bin.
CLTune [40] was the first generic auto-tuning framework with specific support for directly tuning GPU kernels written in the OpenCL programming language. CLTune supports several optimization algorithms, including simulated annealing and particle swarm optimization, but these do not outperform random search [40]. Kernel Tuning Toolkit (KTT) [7] is developed specifically to support online auto-tuning and pipeline tuning, which allows for the exploration of combinations of tunable parameters over multiple kernels. An interesting feature of KTT is the support to keep track of hardware performance counters, such as L2 cache utilization, during benchmarking, which can also be used in advanced search strategies [6]. Auto-Tuning Framework (ATF) [49] implements an innovative way to generate auto-tuning search spaces, for efficient storage and fast exploration of constrained search spaces, but does not focus on introducing new optimization algorithms.
In earlier work, we have introduced Kernel Tuner [61], a generic auto-tuning framework specifically designed to be an easy-to-use and easy to extend tool for researching autotuning optimization algorithms. Kernel Tuner is a state-of-theart framework that implements the largest range of search optimization strategies of all generic auto-tuning frameworks, and was the first generic framework to implement multiple search strategies that consistently outperformed random search [61].

B. Analyzing auto-tuning search spaces
In this paper, we do not only compare the performance of different optimization algorithms on the GPU auto-tuning problem, but we also investigate the properties of the search spaces to understand why certain optimization algorithms outperform others, and gain insight into the difficulty of the optimization problem.
Ryoo et al. [50] were one of the first to study the properties of optimization spaces for GPU applications. They defined two performance metrics to model the efficiency and utilization of a CUDA kernel and used these to find kernel configurations on the pareto curve that maximizes the two metrics. A downside of this approach is that the performance metrics have to be constructed for each kernel individually and require manual inspection and counting instructions in assembly code. Lim et al. [25] also look into search space properties to preselect certain parameter values using static code analysis in order to try to limit exploration of the search space by an auto-tuner.
In [42], Ochoa et al introduce the concept of Local Optima Networks (LONs). When constructing LONs the search space is partitioned into basins of attraction, i.e., sets of points where a local search algorithm will terminate in the same local minimum. A LON is a graph with the local optima as vertices, and a directed edge between nodes if a local search step transforms a solution from one basin of attraction to another. We build upon this idea to define fitness flow graphs (FFGs), which as opposed to LONs contain all the points in the search space. Due to the large number of points with "failure fitnesses" (when a configuration fails to compile) in GPU kernel spaces, defining the basin of attraction is difficult. Instead, we simplify the ideas behind LONs to the entire search space, and quantify how likely local search algorithms terminate in good local optima. To do so, we look at PageRank centrality of local optima. In [15] the idea of using PageRank centrality for LONs as predictor of performance for localseach based heuristics was proposed, and in [14] the PageRank was used to rank space difficulty. We extend this idea to FFGs to determine GPU tuning difficulty.

III. METHOD: OPTIMIZATION PROBLEM
In this section, we define the performance optimization of GPU kernels as a mathematical optimization problem, and we present the optimization algorithms that are part of our experiments.
A. GPU kernels GPU kernels are executed by millions of threads in parallel to perform data-parallel computations on the GPU. However, the compute performance of a GPU kernel depends on how the software has been optimized for the hardware.
There are various different design choices that have an impact on the performance of GPU kernels, and this impact is challenging to accurately predict. For example, the way that a computation is parallelized and mapped on the thread blocks and individual threads affects the utilization of the GPU cores. Other design choices include what data types and data layouts to use in the various memory spaces available to GPU applications. There may also be entirely different algorithms to choose from to implement certain parts of the computation.
Other tunable parameters are introduced through code optimizations that can be enabled and may in turn introduce new parameters, such as tiling factors, vector data types, or partial loop unrolling factors. GPU kernels also have a number inherent parameters in terms of the number of thread blocks and the number of threads per block that are used to execute the kernel. The multitude of implementation choices for GPU kernels result in sizeable, non-convex, and discontinuous kernel design spaces.
To automate the kernel design space exploration process, GPU code can be parameterized, either using a kernel template or a code generator. An auto-tuner can take such a kernel template or code generator and empirically benchmark different kernel configurations, until it has found an efficient implementation. The search performed by the auto-tuner can be treated as a mathematical optimization problem of the form: where f (x) is the performance metric to be minimized, for kernel configuration x for a combination of kernel, GPU device, and input settings. In this work the performance metric to be minimized will be the runtime of the kernel.

B. GPU kernel search spaces
The search space of possible settings for a GPU kernel can be characterized as a finite subset X ⊂ Z n for n different parameters. Specifically, for each dimension 1 ≤ i ≤ n, every entry x i of a point x ∈ X takes values from a finite set S i ⊂ Z. For example, the block dimension might allow values in {16, 32, 64, 128}. The total search space is the Cartesian product of these finite sets The local structure of a search space depends on the definition of neighbouring points. A common definition of the neighbours of a point are the points which differ only in one dimension, and are equal for all other dimensions. Mathematically, according to this definition the set of neighbours N (x) of a point x is Here, we will consider a more restrictive type of neighbourhood concept where we place the additional requirement that the parameter that differs from x i should have a value adjacent to x i in the list S i . For example, if the block dimension is allowed to be [16,32,64,128], then neighbours of x i = 64 would be 32 and 128, and the neighbour of 128 would only be 64. We consider this restriction because this definition gives information on whether closely related parameter values are related in performance.
Points of special interest in the search space are local minima. A point is a local minimum if all neighbouring points have a worse fitness. In other words, there are no improvements to be found in the local neighbourhood. Algorithms that scan local neighbourhoods can get stuck in local minima as there are no close points with better fitness.

C. Black-box optimization algorithms
As the search space is typically too large to iterate over all feasible points, a range of more sophisticated optimization algorithms are used in practice. In this section we describe the optimization algorithms that are considered in the experiments. As a categorization of these algorithms, we distinguish continuous, and discrete algorithms, and algorithms that learn about the stochasticity of the problem.
1) Optimization -discrete algorithms: Since tuning GPU kernels involves choosing the best from a finite set of possibilities, it makes sense to consider discrete optimization algorithms. These algorithms are considered in this work: • Random sampling randomly generates solutions and records the highest scoring one. This strategy serves as a baseline comparison to determine if optimization algorithms offer significant benefits.
We consider several local search (or hill climb) algorithms which iteratively check for a neighbouring solution of lower fitness to visit, until a local minimum is reached. For all local search algorithms we distinguish between best-improvement search, where we move to the best neighbour next, and firstimprovement, where we examine neighbours in a random order and move to one when we encounter an improvement.
Local search algorithms can vary the neighbourhood function that they use to generate new candidates. The algorithms can use both the version outlined in equation 2 (called Hamming), and the more restrictive neighbourhood definition in section III-B (called adjacent). First-improvement variants can decide whether they continue checking the remaining variables first after finding an improvement, or if they restart the search (hyperparameter restart search).
• Multi-start local search (MLS) repeatedly generates random starting solutions, and hill climbs them until a local minimum is reached. Hyperparameters: neighbourhood, restart. • Iterative local search (ILS) [31] is similar to MLS but inherits part of the original local minimum when generating a new starting solution. After reaching a minimum, ILS performs several random permutations to generate a new starting solution. This perturbation size is a tunable parameter. In addition, a tunable exit after no improve hyperparameter randomly restarts if no improvement is found after a that many iterations, which helps to escape basins of attraction for small perturbation sizes. Hyperparameters: perturbation size, exit after no improve, neighbourhood, restart. • Tabu search [9] maintains a queue of previously visited solutions which the algorithm is not allowed to visit. The tunable hyperparameter tabu size defines the queue size, and ensures that the new solution has not been visited for tabu size iterations. Tabu search always picks a new solution, whether it is an improvement or not. Hyperparameters: tabu size, neighbourhood. • Simulated annealing (SA) [21] maintains a temperature parameter that, together with the fitness values, determines the probability that we move to a (potentially worse) neighbouring solution. The temperature parameter is decreased each iteration to mimic the behaviour of cooling processes in material physics. A tunable exploration parameter determines the size of the mutation of the current solution that is performed at each iteration. A hill climber subsequently optimizes the new solution, which is accepted with a certain probability. Hyperparameters: exploration, hill climber, neighbourhood We also consider two discrete population-based algorithms, which require a population size parameter to determine the number of solutions they maintain. Population-based methods iteratively create a next generation of solutions by mixing solutions from the previous generation. A reproduction operator creates new initial solutions from existing ones, e.g., with two-point crossover a section of the solution vector is swapped between two solutions. After new solutions have been created, and their fitnesses determined, a selection mechanism determines which solutions are kept for this generation. For example, in tournament selection a number of randomly picked solutions compete for a spot in the next generation.
• Genetic local search (GLS) [17] (or memetic algorithm) is a population-based method where every solution in a generation is subsequently hill climbed. The initial population is made up of randomly generated solutions, which after hill climbing are all local minima. Next, a number of children is created by reproduction, and a new initial starting population is selected from the batch. These solutions are subsequently hill climbed and the procedure is repeated. Hyperparameters: hill climber, population size, reproduction, selection. • Genetic algorithm (GA) [34] is similar to GLS, but instead of hill climbing each solution, it performs a single mutation only, e.g., permuting a single parameter. Instead of a hill climbing algorithm, GA has a tunable hyperparamter mutation that determines the fraction of variables of a solution that are mutated each generation. Hyperparameters: mutation, population size, reproduction, selection.
2) Optimization -continuous algorithms: As an alternative to discrete algorithms, we can consider continuous optimization algorithms which operate on real-valued solutions. In order to apply algorithms which assume continuous variables to a discrete problem such as GPU kernel tuning, we need to define a mapping between a real-valued vector, and the discrete values in the search space. Suppose the search space allows values x 1 , x 2 , . . . , x n for a particular variable, then a continuous variable y ∈ [0, 1] gets mapped to the closest grid  pointȳ: Effectively, this ensures that all possible discrete values are equally spaced across the interval [0, 1], and the continuous variable is mapped to the closest one. The continuous optimization algorithms operate on real-valued vectors with dimensions equal to the number of parameters that are to be optimized. Each entry is bounded to the unit interval. The mapping ensures that points close together in realvalued space can get mapped to the same point in the GPU tuning space. While this can negatively impact the performance of continuous algorithms, it does not automatically lead to poor performance, as illustrated by the strong performance of continuous algorithms in Kernel Tuner [61]. Furthermore, this mapping allows us to explore a new class of algorithms. Here, we consider two local search algorithms.
• Basin hopping [64] is a global stepping algorithm that chooses new starting positions for local minimization. It requires the local minimizer method and a temperature parameter to be chosen. The temperature parameter determines the accept-reject criterion. Currently supported minimization methods are the nonlinear conjugate gradient (CG) [38], simplex (Nelder-Mead) [37], conjugate direction (Powell) [46], L-BFGS-B [4], Constrained Optimization BY Linear Approximation (COBYLA) [47], and Sequential Least Squares Programming (SLSQP) [22] methods. Hyperparameters: minimizer method, temperature. • Dual annealing [60] is an extension of generalized simulated annealing, paired with a local minimization method. It combines global and local search procedures, and it requires users to choose a local minimization method. Hyperparameters: minimizer method.
Lastly, we consider two population-based algorithms.
• Particle swarm optimization (PSO) [19] initializes a number of particles at random in the search space. Each iteration, these particles update their position and velocity. Particles transmit information to a certain number of neighbours, thereby influencing the movement of the other particles.
Hyperparameters: #particles, neighbours evaluated. • Differential evolution [56] is similar to a genetic algorithm, but mixing strategies are based on real-valued solutions. Typically they involve mixing the best solution with a random candidate, and accepting the result with a certain probability. Hyperparameters: mixing method, population size, mutation size, recombination probability.
3) Optimization -tuning algorithms for stochastic optimization: In this survey we consider two state-of-the-art parameter tuning algorithms; • Sequential Model Algorithm Configuration (SMAC) [27] is a random forest-based Bayesian optimization method that is designed for optimization of stochastic problems. However, it can also be used to optimize deterministic problems. SMAC requires the model type of its Bayesian optimizer to be chosen. The gp-mcmc model was significantly slower and worse than gp in the preliminary experiments. We therefore use the gp model type in this work. The acquisition function of the BO is another tunable hyperparameter. Hyperparameters: acquisition function. • Iterated racing (irace) [30] is a statistical approach for selecting the best configuration out of a set of candidates for stochastic optimization problems. After consulting the authors [30], we set firstTest and nbConfigurations as tunable hyperparameters for irace. Hyperparameters: firstTest, nbConfigurations.

IV. IMPLEMENTATION
In this section we comment on certain implementation details for the software developed for this work. The algorithms and analysis tools are implemented in the BlooPy Python package, and the tuning of the GPU kernels is performed by Kernel Tuner.

A. BlooPy and SOTA packages
The algorithms evaluated in this work are implemented in the discrete optimization package BlooPy (BLackbOx Optimization Python) [  functions for converting discrete solutions that are encoded as lists or arrays to bitstrings. Similarly, continuous solutions are mapped to discrete solution vectors using the mapping outlined in section III-C2. BlooPy requires the search space to be finite. Using the bitstring encodings, BlooPy's algorithms can make use of the computationally efficient Python module bitarray which implements fast low-level bitstrings in C. In addition, the algorithms are automatically applicable to benchmark bitstring-based optimization problems such as randomized N k-landscapes [67]. Optimization algorithms in BlooPy maintain a cache of previously visited solutions. This means that a solution that has been visited before does not count as a function evaluation, and instead the cached value is returned. In addition to a variety of optimization algorithms, BlooPy implements several search space analysis tools. For example, it implements functions to determine the type of points in the search space, e.g., local minima and saddle points. Furthermore, BlooPy implements functions to compute the fitness flow graphs outlined in section VII-B. BlooPy can be installed from the GitHub source repository [51], or by package manager. To perform experiments with SMAC we used the Python package [26], and for irace we used the R package [32].

B. Kernel Tuner
Kernel Tuner [61] implements a wide range of optimization algorithms, and builds on top of various backends (e.g. PyOpenCL, PyCUDA, Cupy, GCC) that take care of the compilation process.
Kernel Tuner runs Python code, provided by the user, which calls the tuner function. In addition, the user needs to provide a code generator or parameterized template for the kernel they wish to optimize. An optimization algorithm then selects different kernel configurations for benchmarking.

V. EXPERIMENTAL SETUP
To analyze the structure of different kernel spaces, and find the optimization algorithm best suited to finding strong kernel settings, we run experiments on 9 different GPUs, for 3 realworld applicable kernel programs (26 kernels in total).
• Convolution [62] operations are an essential tool in image processing, and are often used for tasks such as edge detection, blurring, or sharpening. They also feature prominently in deep learning methods for image processing as they form the backbone of the convolutional neural network (CNN). • GEMM (Generalized dense matrix-matrix multiplication) [39] is one of the most widely-used kernels across many application domains, including neural networks.
Here we perform the calculation C = αA · B + βC for 4096 × 4096 matrices A, B, C, and constants α and β. • PnPoly (Point-in-Polygon) kernel is used by Goncalves et al. [10] as part of a geospatial database management system to, for example, return all objects within the outline of a specific area. Some statistics on the kernel spaces is given in Table II. For convolution and GEMM the majority of the points in the kernel space fail to compile, 68% and 78% respectively. In the case of a failed compilation we attribute a "fail" fitness of 10 10 to this point. The exact kernel spaces can be found in the table of tunable parameters (Table V) in Appendix IX-A.
We selected these kernels programs since they are tunable common subroutines in real-world applications, but also have a compact parameter space which can be fully explored, given ample computation time. Note that this is not feasible for many other kernels used in practice (see section II-A). We have generated cache files of the entire search space for each kernel by brute-force calculation. This allows us to know the optimal settings for each problem, and therefore score solutions returned by algorithms. It also allows us to develop analysis metrics on the entire search space, which could at a later stage be adapted to work when sampling only small parts of the space. We supply our cache files for benchmarking optimization algorithms [52], similar to other computationally expensive applications such as neural architecture search [5].
The 9 GPUs that are used for testing are given in Table I. The convolution kernel is implemented in CUDA, GEMM in OpenCL, and PnPoly is a heterogeneous kernel that runs partly on the CPU, and partly on the GPU using CUDA. The PnPoly kernel uses CUDA-specific features that are not available on AMD GPUs. We have used CUDA Version 11.2, OpenCL 1.2, Python 3.8.5, PyCUDA v2021.1, PyOpenCL v2020.3.1, and BlooPy version 0.4.2.
Experimental setup: The GPU kernels are tuned with respect to runtime (ms). The runtime of a GPU kernel is stochastic, and can vary slightly per execution. Kernel Tuner automatically benchmarks a given configuration 32 times to acquire a mean runtime per configuration. In most cases, the compilation time for a given kernel configuration significantly exceeds the time needed to benchmark 32 runs. Therefore, most of our experiments are performed in a deterministic setting where the fitness of a configuration is the mean runtime. However, we also perform a stochastic experiment where a single kernel runtime is returned for every evaluation. This means that the fitness for the same point in the search space can vary, and algorithms that learn stochastic information, such as irace and SMAC, can potentially benefit.
After discussion with the authors [30] we decided to benchmark irace only for the stochastic experiment. This was decided as it was deemed inappropriate for the deterministic setting since the point of using irace is to dynamically handle stochasticity in expensive problems.
The algorithms are evaluated based on the fraction of the optimal runtime they can find within a limited budget of evaluations. For each algorithm, we run experiments with a  In this section, we first discuss how to initialize each algorithm with favourable hyperparameters. Next, we discuss which algorithms are best suited to tuning GPU kernels.

A. Setting hyperparameters
In order to compare the optimization algorithms fairly for the GPU tuning problem, we need to choose sensible hyperparameters. Which hyperparameters can be varied per algorithm is outlined in italics in section III-C. We test different combinations of hyperparameters on the P100, RTX 2070 Super, and GTX 1080Ti, for all three kernels (9 out of 26 kernels). This way we select reasonable hyperparameters across various architectures and kernels, which users can use as defaults for new GPU tuning problems. We performed a bruteforce search over all combinations of parameter values, and ran each set 20 times for each algorithm. For PSO we kept the w, c 1 , c 2 , p parameters constant as these appeared to have little effect on algorithm performance. Note that the time required to tune hyperparameters varies greatly between algorithms due to this exhaustive search. The hyperparameters chosen are in Table III.
To choose hyperparameters, we first group settings which perform similarly statistically, and attempt to find one set of hyperparameters that performs well across all 9 kernels. For a budget p, let f p,best be the lowest average fitness achieved for a set of hyperparameters, and σ p,best the standard deviation. We perform the following selection approach: 1) For every kernel, we create a set of hyperparameter settings whose average found fitness is within k · σ p,best of f p,best . 2) For each budget, intersect the acceptable settings for convolution, GEMM, and PnPoly, across the 3 GPUs.

3) For each budget, if this intersection is non-empty, reduce
k and repeat. If the intersection is empty, increase k and repeat. Repeat until only one set of hyperparameters remains in the intersection.

B. Kernel tuning algorithm comparison
To quantify which kernel tuning algorithms perform best for certain budgets, we can check whether an algorithm provided statistically significantly better results than others for a certain experiment. To do so, we use a two-sample independent t-test with α = 0.05. For each GPU, kernel, and budget combination, we perform the t-test to see if an algorithm A performed significantly better than algorithm B. We subsequently combine the total number of "wins" for algorithm A across all GPUs (excluding those that were used for tuning).
We split the competitions into low-range (200 evaluations or fewer), and medium-range budgets. The competition tables at different splits can be found in the Appendix IX-B. The full plots per GPU and algorithm can be found in Appendix IX-C. The results of these inter-algorithm competitions are given in the heatmaps displayed in Figure 1. The competition heatmaps display how often the column algorithm found a statistically better solution than the row algorithm in that budget range.   beating basin hopping 5 times, and basin hopping beating dual annealing 6 times. For PnPoly dual annealing was again best, followed by SMAC. We show the total number of wins and losses across all kernels and GPUs in Table IV. Here we see that for low budgets dual annealing has significantly more wins, and fewer losses than all other algorithms. Interestingly, SMAC was second for PnPoly, in the best half of algorithms for convolution, but did not even beat random sampling for GEMM. We hypothesize this is either due to the number of variables to optimize for each of the three kernels (see Table II), or the increasing fraction of fail fitnesses for these kernels. We hypothesize that the Bayesian optimizer could not fit a proper surrogate for GEMM with a low budget, and many failed compilations.
Medium budget: For more than 200 function evaluations, FirstMLS and GLS performed best for convolution, followed by FirstILS. For GEMM, FirstILS and simulated annealing performed best, followed by GLS. For PnPoly, FirstMLS is the strongest algorithm, followed by simulated annealing and GLS. As can be seen from Table IV, FirstILS and simulated annealing have the most number of total wins, and simulated annealing and genetic local search have the least number of losses.
Additional remarks: In general, best-improvement local search algorithms performed significantly worse than the firstimprovement variants. In fact, for the low range, they proved statistically worse than random sampling for PnPoly, and in general have fewer wins and more losses. This can be explained due to the fact that exploring all the neighbours before taking a step costs many evaluations, and leads to only exploring a single neighbourhood for low budgets. For the population-based methods, GLS is the best performing algorithm for medium budgets, but does worse than differential evolution and GA for low budgets. PSO performed significantly worse.
2) Kernel tuning: stochastic fitness: For the stochastic experiments the algorithms have to minimize the runtime of a GPU kernel where the runtime is a random draw from the 32 timings. In addition to running SMAC and irace, we also run FirstILS, GA and dual annealing which did well for deterministic fitnesses. We remark that irace throws an error if the budget is too small with respect to the number of variables, and therefore starts at a budget of 200 for convolution and PnPoly, and 400 for GEMM.
The experimental results are shown in Figure 2. Here we aggregate the results per kernel for all GPU models by showing the mean fraction of optimum (and 95% confidence interval) for a given max budget. We see that GA and dual annealing are best for low budgets in the stochastic experiments. FirstILS does well for budgets ≥ 100. Irace is the best method for GEMM with budgets ≥ 800, but for convolution and pnpoly irace is not as good as GA, dual annealing, and FirstILS.
SMAC consistently achieves a lower fraction of optimality than the competing algorithms across kernels and budgets. Again, we hypothesize that this is because of the high number of fail fitnesses in the search spaces (see Table II). This makes it hard for the Bayesian optimizer to fit a meaningful surrogate.
Stochastic or deterministic: Overall, we notice that higher budgets are necessary to find good solutions for the stochastic experiments than in the deterministic case. This leads to a higher overall tuning time. We therefore recommend to treat GPU kernel tuning as a deterministic optimization problem, with the mean runtime as fitness. The added stochastic information does not appear to allow SMAC or irace to consistently outperform conventional black-box algorithms. This could be because the runtime does not vary much; the average (normalized) runtime and standard deviation is 1.000 ± 0.011. Second, the high number of failure configurations could confuse models that try to learn stochastic information.

VII. QUANTIFYING GPU TUNING DIFFICULTY
In this section we want to gain insight into the difficulty of the GPU kernel tuning optimization problem, and quantify kernel spaces according to tuning difficulty. When attempting to understand why certain GPU kernel spaces appear difficult to optimize we found that relatively simple metrics do not coincide with our experimental results. We outline discrepancies between an intuitive simple metric and our experimental results, and introduce a novel refined approach that does correlate with our results.
A. Naive metric: fraction of optimal fitness of local minima Method: As an example of a simple metric that intuitively could explain the results, we consider the fraction of optimal fitness of local minima. For a minimum x i , and optimal fitness f opt , we can consider the fraction of optimal fitness of the minimum f opt /f (x i ). In this case, we divide the global minimal runtime by the runtime of the minima.
Results: In Figure 3 we show a scatter plot of the fraction of optimal fitness for the local minima for the convolution kernel (per GPU model). According to this distribution, the V100 and A100 GPUs have the closest to optimal median fitness for local minima. This means that an algorithm that randomly explores local minima with equal probability will obtain the closest to optimal runtime for these kernels.
Analysis: To empirically check how difficult the GPU kernels are to tune, we can plot the fraction of optimal fitness that optimization algorithms managed to achieve for certain budgets. If f j is the lowest fitness found for a single run for some budget p, a point in the plot is the average over 50 runs computed asf p := (1/50) ·   and 5 we plotf p for dual annealing and FirstILS. We chose dual annealing and FirstILS as they represent the strongest algorithms for low, and medium budgets respectively. We see that for convolution on the A100 GPU both algorithms returned solutions which were furthest away from the optimum, while for the V100 both optimizers return close to optimal solutions for few function evaluations. These observations are opposite to what would be expected on the basis of Figure 3. Hence, the distribution of fitness of the local optima does not properly explain tuning difficulty.

B. Refined approach: Fitness flow graphs and PageRank
Method: A more refined metric to quantify GPU tuning difficulty may be to compute how likely local search algorithms terminate in local minima. For this purpose, we introduce the fitness flow graph (FFG), which contains all points in the search space, and creates a directed edge to a neighbouring point if the neighbour has lower fitness. This means that a random walk across the FFG mimics the behaviour of a randomized first-improvement local search algorithm. The expected proportion of arrivals of each minimum then gives a metric for weighting reachability of each minimum. We show two example FFGs in Figure 6.
To compute the likelihood of arrival per local minima, we compute the PageRank node centrality, which was originally used to determine the relevance of a webpage [3], [43]. Let A G be the adjacency matrix of a directed graph G, rescaled such that each column adds up to 1. Essentially, this means that for every node, the column is a probability vector of visiting  For an FFG, this means that the PageRank value of a local minimum is the probability of arriving in that minimum after a long random walk through the graph.
As a measure of difficulty we consider how likely a certain subset of "suitably good" local minima are to be visited by a local search algorithm relative to the rest. Suppose that f opt is the optimal fitness, and let L(X) be the set of local minima of X. Given a proportion p, we take the set of nodes L p (X) consisting of local minima with fitness less than (1 + p)f opt (for minimization problems, otherwise (1 − p)f opt ). For a centrality function c G , we define the p-proportion of centrality Results: The proportion of centrality for strong local minima for each FFG is shown in Figure 7. We calculate the proportion of centrality for different acceptance percentages with respect to the global minimum of p = 0, 1, 2, . . . , 15%.
Revisiting the A100 and V100 convolution kernel comparison, we see that the proportion of centrality matches the experimental observations for dual annealing and FirstILS. Figure 7 shows that the NVidia V100 has the most central local minima, whereas the A100 has the least central local minima. For the GEMM and PnPoly kernels, Figures 4 and 5 align with the expectations based on the proportion of centrality. For example, the group of PnPoly kernels with lowest proportion of centrality (P100, GTX Titan X, K20, GTX 1080Ti) are indeed the hardest to tune for both algorithms.
One exception is the K20 GEMM kernel, where proportion of centrality does not entirely reflect the perceived difficulty for dual annealing. This suggests that the proportion of centrality may correlate with GPU tuning difficulty better for certain optimization algorithms. This is to be expected since the PageRank centrality on the FFG in expectation mimics the performance of randomized first-improvement local search. Algorithms that are substantially different than firstimprovement local search will therefore also correlate less with the expected difficulty on the basis of proportion of centrality.
Analysis: Overall, the experimental results suggest that the proportion of centrality is a suitable metric for estimating tuning difficulty for GPU kernels. By using FFGs and the PageRank algorithm, we are able to observe kernel differences that were otherwise unknown. For example, both the A100 and V100 convolution kernels have few outlier minima with a close to optimal fitness. In fact, the existence of only a few kernel configurations that lead to large increases in performance is a general property of certain GPU kernels [61]. Crucially however, the likelihood of local search algorithms arriving in such minima differs greatly between the A100 and V100. The proportion of centrality of an FFG gives us a tool to quantify this likelihood. However, further research is necessary to quantitatively determine how well our proposed metric correlates with GPU tuning difficulty.
As a final remark on kernel differences, the experimental results shows that the difficulty of tuning a particular kernel can greatly differ from one GPU to the next, and that these changes do not appear to be correlated with release time of the models. The A100 is the most recent GPU in our set, while the K20 is the oldest. For GEMM and PnPoly, we can say that it has become easier to tune these kernels with more recent GPUs, but the convolution kernel has become more difficult to tune, except on the V100.

VIII. CONCLUSION
In this paper, we have investigated which optimization algorithms produce the fastest GPU kernel configurations across different tuning-time ranges. To do so, we analyzed 26 GPU kernel spaces for 9 GPUs. We computed sets of optimal hyperparameters for GPU tuning for each optimization algorithm. From among the tested algorithms in this set of experiments, we conclude that dual annealing performs best as GPU kernel tuner when a limited amount of function evaluations is desirable. When more evaluations are possible, firstimprovement local searchers such as FirstILS proved the best GPU kernel tuners. Using these algorithms, we are convinced that GPU programmers can reliably auto-tune GPU kernels to close to optimal runtime while requiring relatively few re-compilations of the code. Furthermore, we conclude that   treating GPU tuning as a deterministic optimization problem is preferred over treating the runtime as a stochastic variable.
We showed that the basic metric of fraction of optimality of local minima is not suitable for explaining the results observed in the experimental benchmarks. To make steps towards a metric for tuning difficulty, we introduced the concept of fitness flow graphs, and proportion of centrality. Our results suggest that the proportion of centrality can be used to quantify tuning difficulty. For future work, in cases where exhaustive exploration is infeasible, perhaps a procedure to dynamically update the proportion of centrality of FFGs can be used. Such dynamic estimates of tuning difficulty could be used for automatic algorithm selection within frameworks such as Kernel Tuner. Furthermore, the pagerank centrality of strong local minima within FFGs can be used to investigate why certain minima are unlikely to be visited, for example because neighbouring configurations fail to compile. Lastly, in this work we fully computed 26 kernel spaces, and made these publicly available. We aim to extend this to a benchmark dataset for evolutionary computation algorithms.  IX. APPENDIX A. Appendix: Tunable parameters per GPU kernel In Table V we show the tunable parameters per kernel, and the values each parameter could take. For the convolution kernel, the MI50 GPU (the only AMD model) required a different problem setup due to hardware constraints.

B. Appendix: Alternative splits for competition heatmaps
In Figures 8 and 9 we show the algorithm competition heatmaps such as in Figure 1 but when split at 100 and 400 budgets instead of 200.