# Agent-based Collaborative Random Search for Hyper-parameter Tuning and
Global Function Optimization

Ahmad Esmaeili Corresponding Author Zahra Ghorrati  Eric T. Matson

( 1Department of Computer and Information Technology, Purdue University, West
Lafayette, IN 47907
{aesmaei, zghorrat<EMAIL_ADDRESS>
)

###### Abstract

Hyper-parameter optimization is one of the most tedious yet crucial steps in
training machine learning models. There are numerous methods for this vital
model-building stage, ranging from domain-specific manual tuning guidelines
suggested by the oracles to the utilization of general-purpose black-box
optimization techniques. This paper proposes an agent-based collaborative
technique for finding near-optimal values for any arbitrary set of hyper-
parameters (or decision variables) in a machine learning model (or general
function optimization problem). The developed method forms a hierarchical
agent-based architecture for the distribution of the searching operations at
different dimensions and employs a cooperative searching procedure based on an
adaptive width-based random sampling technique to locate the optima. The
behavior of the presented model, specifically against the changes in its
design parameters, is investigated in both machine learning and global
function optimization applications, and its performance is compared with that
of two randomized tuning strategies that are commonly used in practice.
According to the empirical results, the proposed model outperformed the
compared methods in the experimented classification, regression, and multi-
dimensional function optimization tasks, notably in a higher number of
dimensions and in the presence of limited on-device computational resources.

Keywords: Multi-Agent Systems; Distributed Machine Learning; Hyper-Parameter
Tuning; Agent-based Optimization; Random Search

##  1 Introduction

Almost all Machine Learning (ML) algorithms comprise a set of hyper-parameters
that control their learning process and the quality of their resulting models.
The number of hidden units, the learning rate, the mini-batch sizes, etc. in
neural networks; the kernel parameters and regularization penalty amount in
support vector machines; and maximum depth, samples split criteria, and the
number of used features in decision trees are few common hyper-parameter
examples that need to be configured for the corresponding learning algorithms.
Assuming a specific ML algorithm and a dataset, one can build countless number
of models, each with potentially different performance and/or learning speeds,
by assigning different values to the algorithm’s hyper-parameters. While they
provide ultimate flexibility in using ML algorithms in different scenarios,
they also account for most failures and tedious development procedures.
Unsurprisingly, there are numerous studies and practices in the machine
learning community devoted to the optimization of hyper-parameter. The most
straightforward yet difficult approach utilizes expert knowledge to identify
potentially better candidates in hyper-parameter search spaces to evaluate and
use. The availability of expert knowledge and generating reproducible results
are among the primary limitation of such manual searching technique [1],
particularly due to the fact that using any learning algorithm on different
datasets likely requires different sets of hyper-parameter values [2].

Formally speaking, let $\Lambda=\\{\bm{\lambda}\\}$ denote the set of all
possible hyper-parameter value vectors and
$\mathcal{X}=\\{\mathcal{X}^{(train)},\mathcal{X}^{(valid)}\\}$ be the dataset
split into training and validation sets. The learning algorithm with hyper-
parameter values vector $\bm{\lambda}$ is a function that maps training
dataset $\mathcal{X}^{(train)}$ to model $\mathcal{M}$, i.e.
$\mathcal{M}=\mathcal{A}_{\bm{\lambda}}(\mathcal{X}^{(train)})$, and the
hyper-parameter optimization problem can be formally written as [1]:

$\bm{\lambda}^{(*)}=\operatorname*{arg\,min}_{\bm{\lambda}\in\Lambda}\mathbb{E}_{x\sim\mathcal{G}_{x}}\left[\mathcal{L}\left(x;\mathcal{A}_{\bm{\lambda}}(\mathcal{X}^{(train)})\right)\right]$
(1)

where $\mathcal{G}_{x}$ and $\mathcal{L}(x;\mathcal{M})$ are respectively the
grand truth distribution and the expected loss of applying learning model
$\mathcal{M}$ over i.i.d. samples $x$; and
$\mathbb{E}_{x\sim\mathcal{G}_{x}}\left[\mathcal{L}\left(x;\mathcal{A}_{\bm{\lambda}}(\mathcal{X}^{(train)})\right)\right]$
gives the generalization error for algorithm $\mathcal{A}_{\bm{\lambda}}$. To
cope with the inaccessibility of the grand truth in real-world problems, the
generalization error is commonly estimated using the cross-validation
technique [3], leading to the following approximation of the above-mentioned
optimization problem:

$\bm{\lambda}^{(*)}\approx\operatorname*{arg\,min}_{\bm{\lambda}\in\Lambda}\operatorname*{mean}_{x\in\mathcal{X}^{(valid)}}\mathcal{L}\left(x;\mathcal{A}_{\bm{\lambda}}(\mathcal{X}^{(train)})\right)\equiv\operatorname*{arg\,min}_{\bm{\lambda}\in\Lambda}\Psi(\bm{\lambda})$
(2)

where $\Psi(\bm{\lambda})$ is called hyper-parameter response function [1].

Putting the manual tuning approaches aside, there is a wide range of
techniques that use global optimization methods to address the ML hyper-
parameter tuning problem. Grid search [4, 5], random search [1], Bayesian
optimization [6, 7, 8], and evolutionary and population-based optimization [9,
10, 11] are some common class of tuning methodologies that are studied and
used extensively by the community. In grid search, for instance, every
combination of a predetermined set of values in each hyper-parameter is
evaluated and the hyper-parameter value vector that minimizes the loss
function is selected. For $k$ number of configurable hyper-parameters, if we
denote the set of candidate values for the $j$-th hyper-parameter
$\lambda_{j}^{(i)}\in\bm{\lambda}^{(i)}$ by $\mathcal{V}_{j}$, grid search
would evaluate $T=\Pi_{j=1}^{k}|\mathcal{V}_{j}|$ number of trials that can
grow exponentially with the increase in the number of configurable hyper-
parameters and the quantity of the candidate values for each dimension. This
issue is referred to as the curse of dimensionality [12] and is the primary
reason for making grid search an uninteresting methodology in large-scale
real-world scenarios. In the standard random search, on the other hand, a set
of $b$ uniformly distributed random points in the hyper-parameter search
space, $\\{\lambda^{(1)},\dots,\lambda^{(b)}\\}\in\Lambda$, are evaluated to
select the best candidate. As the number of evaluations only depends on the
budget value, $b$, random search does not suffer from the curse of
dimensionality, is shown to be more effective than grid search [1], and is
often used as a baseline method. Bayesian optimization, as a global black box
expensive function optimization technique, iteratively fits a surrogate model
to the available observations, $(\lambda^{(i)},\Psi_{\lambda^{(i)}})$ and
then, uses an acquisition function to determine the next hyper-parameter
values to evaluate and use in the next iteration [8, 13]. Unlike grid and
random search methods, in which the searching operations can be easily
parallelized, the Bayesian method is originally sequential though various
distributed versions have been proposed in the literature [14, 15].
Nevertheless, thanks to its sample efficiency and robustness to noisy
evaluations, Bayesian optimization is a popular method in the hyper-parameter
tuning of deep-learning models, particularly when the number of configurable
hyper-parameters is less than 20 [16]. Evolution and population-based global
optimization methods, such as genetic algorithms and swarm-based optimization
techniques, form the other class of common tuning approaches, in which the
hyper-parameters configurations are improved over multiple generations
generated by local and global perturbations [17, 10, 11]. Population-based
methods are embarrassingly parallel [18] and, similar to grid and random
search approaches, the evaluations can be distributed over multiple machines.

Multi-Agent Systems (MAS) and agent-based technologies, when applied to
machine learning and data mining, bring about scalability and autonomy, and
facilitate the decentralization of learning resources and the utilization of
strategic and collaborative learning models [19, 20, 21]. Agnet-based machine
learning and collaborative hyper-parameter tuning are not new to this paper
and have previously been studied in the literature. The research reported in
[22] is among the noteworthy contributions, which proposes a surrogate-based
collaborative tuning technique incorporating the experience achieved from
previous experiments. To put it simply, this model performs simultaneous
configuration of the same hyper-parameters over multiple datasets and employs
the gained information in all subsequent tuning problems. Auto-Tuned Models
(ATM) [23] is a distributed and collaborative system that automates hyper-
parameter tuning and classification model selection procedures. At its core,
ATM utilizes Conditional Parameter Tree (CPT), in which a learning method is
placed at the root and its children are the method’s hyper-parameters, to
represent the hyper-parameter search space. Different tunable subsets of
hyper-parameter nodes in CPT are selected during model selection, and assigned
to a cluster of workers to be configured. Koch et al. [24] introduced Autotune
as a derivative-free hyper-parameter optimization framework. Being composed of
a hybrid and extendable set of solvers, this framework concurrently runs
various searching methods potentially distributed over a set of workers to
evaluate objective functions and provide feedback to the solvers. Autotune
employs an iterative process during which all the points that have already
been evaluated are exchanged with the solvers to generate new sets of points
to evaluate. In learning-based settings, the work reported in [25] used Mutli-
Agent Reinforcement Learning (MARL) to optimize the hyper-parameters of deep
convolutional neural networks (CNN). The suggested model splits the design
space into sub-spaces and devotes each agent to tuning the hyper-parameters of
a single network layer using Q-learning. Parker-Holder et al. [26] presented
the Population-Based Bandit (PB2) algorithm, which efficiently directs the
searching operation of hyper-parameters in Reinforcement Learning using a
probabilistic model. In PB2, a population of agents is trained in parallel,
and their performance is monitored on a regular basis. An underperforming
agent’s network weights are replaced with those of a better-performing agent,
and its hyper-parameters are tuned using Bayesian optimization.

In continuation to our recent generic collaborative optimization model [21],
this paper delves into the design of a multi-level agent-based distributed
random search technique that can be used for both hyper-parameter tuning and
general-purpose black-box function optimization. The proposed methods, at its
core, forms a tree-like structure comprising a set of interacting agents that,
depending on their position in the hierarchy, focus on tuning/optimizing a
single hyper-parameter/decision variable using a biased hyper-cube random
sampling technique, or aggregating the results and facilitating collaborations
based on the gained experience of other agents. The rationales behind choosing
random search as the core tuning strategy of the agents include, but are not
limited to, its intrinsic distributability, acceptable performance in
practice, and not requiring differentiable objective functions. Although the
parent model in [21] does not impose any restrictions on the state and
capabilities of the agents, this paper assumes homogeneity in the sense that
the tuner/optimizer agents use the same mechanism for their assigned job. With
that said, the proposed method is analyzed in terms of its design parameters,
and the empirical results from conducted ML classification and regression
tasks, as well as, various multi-dimensional functions optimization problems
demonstrate that the suggested approach not only outperforms the underlying
random search methodologies under the same deployment conditions but also
provide a better-distributed solution in the presence of limited computational
resources.

The remainder of this paper is organized as follows: section 2 dissects the
proposed agent-based random search method; section 3 presents the details of
used experimental ML and function optimization settings, and discusses the
performance of the proposed model under different scenarios; and finally,
section 4 concludes the paper and provides future work suggestions.

##  2 Methodology

This section dissects the proposed agent-based hyperparameter tuning and
black-box function optimization approach. To help with a clear understanding
of the proposed algorithms, this section begins by providing the preliminaries
and introducing the key concepts and then presents the details of the agent-
based randomized search algorithms accompanied by hands-on examples whenever
needed.

###  2.1 Preliminaries

An agent, as the first-class entity in the proposed approach, might play
different roles depending on its position in the system. As stated before,
this paper is using hierarchical structures to coordinate the agents in the
system and hence, it defines two agent types: (1) internals, which play the
role of result aggregators and collaboration facilitators in connection with
their subordinates; and (2) terminals, which, implementing a single-variable
randomized searching algorithm, are the actual searchers/optimizers positioned
at the bottom-most level of the hierarchy. Assuming $G$ to be the set of all
agents in the system and the root of the hierarchy to be at level 0, this
paper uses $g_{\lambda_{i}}^{l}$ ( $G_{\bm{\lambda_{j}}}^{l}$) to denote the
agent (the set of agents) at level $l$ of the hierarchy that are specialized
in tuning hyper-parameter $\lambda_{i}$ (hyper-parameter set
$\bm{\lambda_{j}}$) respectively, where
$\bm{\lambda_{j}}\subseteq\bm{\lambda}$ and $g_{\lambda_{i}}^{l+1}\in
G_{\bm{\lambda_{j}}}^{l}$ iff. $\lambda_{i}\in\bm{\lambda_{j}}$.

As denoted above, the hyper-parameters that the agents represent determine
their position in the hierarchy. let
$\mathcal{A}_{\bm{\lambda}=\\{\lambda_{2},\lambda_{2},\ldots\lambda_{n}\\}}$
be the ML algorithms for which we intend to tune the hyper-parameters. As the
tuning process might not target the entire hyper-parameter set of the
algorithm, the proposed method divides the set into two objective and fixed
disjoint subsets, which respectively denoted by $\bm{\lambda_{o}}$ and
$\bm{\lambda_{f}}$ refer to the hyper-parameter sets that we intend to tune
and the ones we need to keep fixed. Formally, that is
$\bm{\lambda}=\bm{\lambda_{o}}\cup\bm{\lambda_{f}}$ and
$\bm{\lambda_{o}}\cap\bm{\lambda_{f}}=\emptyset$. The paper further assumes
two types of objective hyper-parameters: (1) primary hyper-parameters, denoted
by $\bm{\hat{\lambda}_{o}}$, which comprise the main targets of the
corresponding tuners (agents); and (2) subsidiary hyper-parameters, denoted by
$\bm{\hat{\lambda}^{\prime}_{o}}$, which include the ones whose values are set
by the other agents to help limit the searching space. These two sets are
complements of each other, i.e.,
$\bm{\hat{\lambda}^{\prime}_{o}}=\bm{\lambda_{o}}-\bm{\hat{\lambda}_{o}}$, and
the skill of an agent is determined by the primary objective set,
$\bm{\hat{\lambda}_{o}}$, it represents. With that said, for all terminal
agents in the hierarchy, we have $|\bm{\hat{\lambda}_{o}}|=1$, where $|\dots|$
denotes the set cardinality.

The agents of a realistic MAS are susceptible to various limitations that are
imposed by their environment and/or computational resources. This paper, due
to its focus on the decentralization of the searching process, foresees two
limitation aspects for the agents: (1) the maximum number of concurrent
connections, denoted by $c$, that an agent can manage; (2) the number of
concurrent processes, called budget and denoted by $b$, that the agent can
execute and handle. In the proposed method, $c>1$ determines the maximum
number of subordinates (children) that an internal agent can have. The budget,
$b\geq 1$, on the other hand, put a restriction on the maximum number of
parallel evaluations that an agent can perform in the step of searching for
the optima.

Communications play a critical role in all MAS, including the agent-based
method proposed in this paper. For all intra-system, i.e. between any two
agents, and inter-system, i.e. between an agent and a user, interactions, the
suggested method uses the following tuple-based structure for the queries:

$\left<\mathcal{A}_{\bm{\lambda}},\\{\bm{\hat{\lambda}_{o}},\bm{\hat{\lambda}^{\prime}_{o}},\bm{\lambda_{f}}\\},\bm{\mathcal{V}},\\{\mathcal{X}^{(train)},\mathcal{X}^{(valid)}\\},\mathcal{L}\right>$
(3)

where, $\bm{\mathcal{V}}=\\{(\lambda_{i},v_{i})\\}_{1\leq i\leq n}$ denotes
the set containing the candidate values for all hyper-parameter, and the
remaining notations are as defined in equation 1. Based on what discussed
before, it is clear that
$|\bm{\lambda_{f}}|\leq|\bm{\mathcal{V}}|\leq\bm{\lambda}$.

###  2.2 Agent-based Randomized Searching Algorithm

The high-level abstract view of the proposed approach is composed of two major
steps: (1) distributedly building the hierarchical MAS; and (2) performing the
collaborative searching process through vertical communications in the
hierarchy. The sections that follow go into greater detail about these two
stages.

####  2.2.1 Distributed Hierarchy Formation

As for the first step, each agent, $t$, divides the primary objective hyper-
parameter set of the query it receives, i.e. $\bm{\hat{\lambda}_{o}}$, into
$c_{t}>1$ number of subsets, for each of which, the system initiates a new
agent to handle. This process continues recursively until there is only one
hyper-parameter in the primary objective set, i.e.
$|\bm{\hat{\lambda}_{o}}|=1$, which is assigned to a terminal agent. Figure 1
provides an example hierarchy resulted from the recursive division of primary
objective set
$\bm{\lambda_{o}}=\\{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4},\lambda_{5},\lambda_{6}\\}$.
For the sake of clarity, we have used the indexes of the hyper-parameters as
the labels of the nodes in the hierarchy, and the green and orange colors are
employed to highlight the members of the
$\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}{\bm{\hat{\lambda}_{o}}}$
and
$\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}{\bm{\hat{\lambda}^{\prime}_{o}}}$
sets respectively. Regarding the maximum number of concurrent connections that
the agents can handle in this example, it is assumed for all agents $c=2$,
except for the rightmost agent in the second level of the hierarchy, for which
$c=3$. It is worth emphasizing that at the beginning of the process, when the
tuning query is received from the user, we have
$\bm{\hat{\lambda}_{o}}=\bm{\lambda_{o}}$ and
$\bm{\hat{\lambda}^{\prime}_{o}=\emptyset}$, which is the reason for the all
green node of the root node in this example.

Figure 1: Hierarchical structure built for
$\bm{\lambda_{o}}=\\{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4},\lambda_{5},\lambda_{6}\\}$,
where the primary and complementary hyper-parameters of each node are
respectively highlighted in green and orange, and the labels are the indexes
of $\lambda_{i}$.

Algorithm 1 presents the details of the process. We have chosen self-
explanatory names for the functions and variables and provided comments
wherever they are required to improve clarity. In this algorithm, the function
PrepareResources in line 1 prepares the data and computational resources for
the newly built/assigned terminal agent. Such resources are used for training,
validation, and tuning process. The function SpawnOrConnect in line 1 creates
a subordinate agent that represents ML algorithm $\mathcal{A}_{\lambda}$ and
expected loss function $\mathcal{L}$. This is done by either creating a new
agent or connecting to an existing idle one if the resources are reused. Two
functions PrepareFeedback and Tune in lines 1 and 1, respectively, are called
when the structure formation process is over and the root agent initiates the
tuning process in the hierarchy. Later, these two functions are discussed in
more detail.

Algorithm 1 Distributed formation of the hierarchical agent-based hyper-
parameter tuning structure.

Function  _Start(_$\left
<\mathcal{A}_{\bm{\lambda}},\\{\bm{\hat{\lambda}_{o}},\bm{\hat{\lambda}^{\prime}_{o}},\bm{\lambda_{f}}\\},\bm{\mathcal{V}},\\{\mathcal{X}^{(train)},\mathcal{X}^{(valid)}\\},\mathcal{L}\right>$_)_:

1 if  _$|\bm{\hat{\lambda}_{o}}|=1$_ then $\triangleright$ agent is terminal

2
$\mathcal{R}\leftarrow$PrepareResources($\left<\\{\bm{\hat{\lambda}_{o}},\bm{\hat{\lambda}^{\prime}_{o}},\bm{\lambda_{f}}\\},\\{\mathcal{X}^{(train)},\mathcal{X}^{(valid)}\\}\right>$)
Inform($Parent,\mathcal{R}$) $\triangleright$ informs the parent agent

3 else $\triangleright$ agent is internal ($|\bm{\hat{\lambda}_{o}}|>1$)

4 $k\leftarrow\min(c_{my},|\bm{\hat{\lambda}_{o}}|)$ $\triangleright$ the
number of children
for  _$i\leftarrow 1$ to $k$_ do

5 $G_{i}\leftarrow$ SpawnOrConnect($\mathcal{A}_{\bm{\lambda}},\mathcal{L}$)
$\bm{\hat{\lambda}_{o_{i}}}\leftarrow$Divide($\bm{\hat{\lambda}_{o}}$,$i$,
$k$) $\triangleright$ the $i^{\text{th}}$ unique devision
$\bm{\hat{\lambda}^{\prime}_{o_{i}}}\leftarrow(\bm{\hat{\lambda}_{o}}-\bm{\hat{\lambda}_{o_{i}}})\cup\bm{\hat{\lambda}^{\prime}_{o}}$
$\mathcal{R}_{i}\leftarrow$Ask($G_{i}$, Start,
$\left<\mathcal{A}_{\bm{\lambda}},\\{\bm{\hat{\lambda}_{o_{i}}},\bm{\hat{\lambda}^{\prime}_{o_{i}}},\bm{\lambda_{f}}\\},\bm{\mathcal{V}},\\{\mathcal{X}^{(train)},\mathcal{X}^{(valid)}\\},\mathcal{L}\right>$)

6 end for

7 $\mathcal{R}\leftarrow$Aggregate($\\{\mathcal{R}_{i}\\}_{i}$)
$\triangleright$ combines children’s answers
if  _Parent $\neq\emptyset$_ then

8 Inform($Parent,\mathcal{R}$)

9 else

10
$\mathcal{F}\leftarrow\textsc{PrepareFeedback}(\mathcal{R},\bm{\mathcal{V}})$
Tune($\mathcal{F}$) $\triangleright$ initiates the tuning process

11 end if

12

13 end if

14

15 end

####  2.2.2 Collaborative tuning Process

The collaborative tuning process is conducted through a series of vertical
communications in the built hierarchy. Initiated by the root agent, as
explained in the previous section, the Tune request is propagated to all of
the agents in the hierarchy. As for the internal agents, the request will be
simply passed down to the subordinates as they arrive. As for the terminal
agents, on the other hand, the request launches the searching process in the
sub-space specified by the parent. The flow of the results will be in an
upward direction with a slightly different mechanism. As soon as a local
optimum is found by a terminal agent, it will be sent up to the parent agent.
Having waited for the results to be collected from all of its subordinates,
the parent aggregates them together and passes the combined result to its own
parent. This process continues until it reaches the root agent where the new
search guidelines are composed for the next search round.

Algorithm 2 presents the details of the iterated collaborative tuning process,
which might be called by both terminal and internal agents. When it is called
by a terminal agent, it initiates the searching operation for the optima of
the hyper-parameter that the agent represents and informs the result to its
parent. Let $g_{\lambda_{j}}^{l}$ be the terminal agent concentrating on
tuning hyper-parameter $\lambda_{j}$. As it can be seen in line 2, the result
of the search will be a single-item set composed of the identifier of the
hyper-parameter, i.e. $\lambda_{j}$, the set $\bm{\mathcal{V}_{j}}^{(*)}$
containing the coordinate of the best candidate agent $g_{\lambda_{j}}^{l}$
has found, and the response function value for that best candidate, i.e.
$\Psi^{(*)}_{j}$. An internal agent running this procedure, on the other hand,
merely passes the tuning request to the subordinates and waits for their
search results (line 2 of the algorithm). Please note that this asking
operation comprises filtering operation on set $\bm{\mathcal{F}}$. That is,
the subordinate will receive a subset
$\bm{\mathcal{F}}_{i}\subset\bm{\mathcal{F}}$ that only includes starting
coordinates for the terminal agents reachable through that agent. Having
collected all of the results from its subordinates, the internal agent
aggregates them by simply joining the result sets and informs its own parent,
in case it is not the root agent. This process is executed recursively until
the aggregated results reach the root agent of the hierarchy. Depending on
whether the stopping criteria of the algorithm are reached, the root prepares
feedback to initiate the next tuning iteration or a report detailing the
results. The collaboration between the agents is conducted implicitly through
the feedback that the root agent provides for each terminal agent based on the
results it has gathered in the previous iteration. As presented in line 1 of
the algorithm, this feedback basically determines the coordinates of the
position where the terminal agents should start their searching operation. It
should be noted that the argmin function in this operation is due to employing
loss function $\mathcal{L}$ as a metric to evaluate the performance of an ML
model. For performance measures in which maximization is preferred, such as in
accuracy, this operation needs to be replaced by argmax accordingly.

Algorithm 2 Iterated collaborative tuning procedure.

Function  _Tune(_$\bm{\mathcal{F}}$_)_:

16 if  _$Children=\emptyset$_ then $\triangleright$ terminal agent agent

17
$\\{(\lambda_{j},\bm{\mathcal{V}}^{(*)}_{j},\Psi^{(*)}_{j})\\}\leftarrow$RunTuningAlgorithm($\bm{\mathcal{F}}=\bm{\mathcal{V}}$)
Inform($Parent,\\{(\lambda_{j},\bm{\mathcal{V}}^{(*)}_{j},\Psi^{(*)}_{j})\\}$)

18 else

19 foreach  _$G_{i}^{l+1}\in{Children}$_ do

20 $\bm{\mathcal{R}}^{(*)}_{i}\leftarrow$Ask($G_{i}^{l+1}$, Tune,
$\bm{\mathcal{F}_{i}}\subset\bm{\mathcal{F}}$) $\triangleright$
$\bm{\mathcal{R}}^{(*)}_{i}=\\{(\lambda_{k},\bm{\mathcal{V}}^{(*)}_{k},\Psi^{(*)}_{k})\\}_{k}$

21 end foreach

22 $\bm{\mathcal{R}}^{(*)}\leftarrow\bigcup\limits_{G_{i}^{l+1}\in
Children}\bm{\mathcal{R}}^{(*)}_{i}$ $\triangleright$ aggregates results
if  _$Parent\neq\emptyset$_ then $\triangleright$ non-root internal agent

23 Inform($Parent,\bm{\mathcal{R}}^{(*)}$)

24 else

25 if  _ShouldStop($StopCriteria$)$\neq\text{True}$_ then

26
$\bm{\mathcal{F}}\leftarrow\left\\{(\lambda_{i},\bm{\mathcal{V}}_{j});j=\operatorname*{arg\,min}\limits_{1\leq
k\leq n}\Psi^{(*)}_{k}\right\\}_{1\leq i\leq n}$ $\triangleright$ prepares
feedback Tune($\bm{\mathcal{F}}$) $\triangleright$ initiates next tuning
iteration

27 else

28 PrepareReport($\bm{\mathcal{R}}^{(*)}$) $\triangleright$ reports final
result

29 end if

30

31 end if

32

33 end if

34

35 end

The details of the tuning function that each terminal agent runs in line 2 of
algorithm 2 to tune a single hyper-parameter are presented in algorithm 3. As
its input, this function receives a coordinate that the agent
$g_{\lambda_{i}}^{l}$ will use as its starting point in the searching process.
The received argument together with $b$ additional coordinates that the agent
generates randomly are stored in the set of candidate $\bm{\mathcal{C}}$.
Accordingly, $\bm{\mathcal{C}}[c]$ and $\bm{\mathcal{C}}[c](\lambda_{i})$
refer to the $c$-th coordinate in the set and the value assigned to hyper-
parameter $\lambda_{i}$ of that coordinate respectively. Moreover, please
recall from section 2.1 that $b$ denotes the evaluation budget of a terminal
agent. The terminal agents in the proposed method employ a slot-based uniform
random sampling to explore the search space. Formally, let
$\bm{\mathcal{E}}=\\{\epsilon_{\lambda_{1}},\epsilon_{\lambda_{2}},\dots,\epsilon_{\lambda_{n}}\\}$
be a set of real values that each agent utilizes for each hyper-parameter to
control the size of slots in any iteration. Similarly, let
$\bm{s}=\\{s_{\lambda_{1}},s_{\lambda_{2}},\dots,s_{\lambda_{n}}\\}$ specify
the coordinate of the position that an agent starts its searching operation in
any iteration. To sample $b$ random values in the domain
$\mathbb{D}_{\lambda_{j}}$ of any arbitrary hyper-parameter $\lambda_{j}$, the
agent will generate one uniform random value in range

$\mathcal{R}=\left[\max(\inf\mathbb{D}_{\lambda_{j}},s_{\lambda_{j}}-\epsilon_{\lambda_{j}}),\min(\sup\mathbb{D}_{\lambda_{j}},s_{\lambda_{j}}+\epsilon_{\lambda_{j}}))\right]$
(4)

and $b-1$ random values in $\mathbb{D}_{\lambda_{j}}-\mathcal{R}$ by splitting
it into $b-1$ slots and choosing one uniform random value in each slot (lines
3 and 3 of the algorithm). The generation of the uniform random values is done
by calling function UniformRand($A_{1},A_{2},A_{3}$). This function divides
range $A_{1}$ into $A_{2}$ equal-sized slots and returns the uniform random
value generated in the $A_{3}$-th slot. As it can be seen in line 3 of the
algorithm, the agent employs the same function to generate one and only one
value per each element in its subsidiary objective hyper-parameter set
$\bm{\hat{\lambda}^{\prime}_{o}}$.

The slot width parameter set $\bm{\mathcal{E}}$ is used to control the
exploration behavior of the agent around the starting coordinates in the
search space. For instance, for any arbitrary hyper-parameter $\lambda_{i}$,
very small values of $\epsilon_{\lambda_{i}}$ emphasize generating candidates
in the close vicinity of the starting position. On the other hand, larger
values of $\epsilon_{\lambda_{i}}$ decrease the chance that the generated
candidate be close to the starting position. In the proposed method, the
agents adjust $\bm{\mathcal{E}}$ adaptively. To put it formally, each agent
starts the tuning process with the pre-specified value set
$\bm{\mathcal{E}}^{(0)}$, and assuming that $\bm{\mathcal{C}}^{(*)}$ denotes
the best candidate that the agent has found in the previous iteration, the
width parameter set $\bm{\mathcal{E}}$ in iteration $i$ is updated as follows:

$\bm{\mathcal{E}}^{(i)}=\begin{cases}\bm{\Delta}\odot\bm{\mathcal{E}}^{(i-1)}&\text{If
}\bm{\mathcal{V}}=\bm{\mathcal{C}}^{(*)}\\\
\bm{\mathcal{E}}^{(i-1)}&\text{otherwise}\end{cases}$ (5)

where
$\Delta=\\{\delta_{\lambda_{1}},\delta_{\lambda_{2}},\dots,\delta_{\lambda_{n}}\\}$
denotes the scaling changes to apply to the width parameters, and $\odot$
denotes the element-wise multiplication operator. As the paper discusses in
section 3, despite the generic definitions provided here for futuristic
extensions, using the same scaling value for all primary hyper-parameters has
led to satisfactory results in our experiments.

Algorithm 3 A terminal agent’s randomized tuning process.

Function
_RunTuningAlgorithm(_$\bm{\mathcal{F}}=\bm{\mathcal{V}}=\\{(\lambda_{m},v_{m})\\}_{1\leq
m\leq n}$_)_:

$\bm{\mathcal{C}}[0]\leftarrow\bm{\mathcal{V}}$
$\mathcal{R}_{\lambda_{i}}\leftarrow\left[\max(\inf\mathbb{D}_{\lambda_{i}},v_{i}-\epsilon_{\lambda_{i}}),\min(\sup\mathbb{D}_{\lambda_{i}},v_{i}+\epsilon_{\lambda_{i}}))\right]$
for  _$c\leftarrow 1$ to $c=b$_ do

36 if  _$c=1$_ then $\triangleright$ the first sample for
$\bm{\hat{\lambda}_{o}}=\\{\lambda_{i}\\}$

37
$\bm{\mathcal{C}}[c](\lambda_{i})\leftarrow$UniformRand($\mathcal{R}_{\lambda_{i}},1,1$)

38 else $\triangleright$ the remaining samples for
$\bm{\hat{\lambda}_{o}}=\\{\lambda_{i}\\}$

39
$\bm{\mathcal{C}}[c](\lambda_{i})\leftarrow$UniformRand($\mathbb{D}_{\lambda_{i}}-\mathcal{R}_{\lambda_{i}},b-1,c-1$)

40 end if

41 forall  _$\lambda_{k}\in\bm{\hat{\lambda}^{\prime}_{o}}$_ do

42
$\mathcal{R}_{\lambda_{k}}\leftarrow\left[\max(\inf\mathbb{D}_{\lambda_{k}},v_{k}-\epsilon_{\lambda_{k}}),\min(\sup\mathbb{D}_{\lambda_{k}},v_{k}+\epsilon_{\lambda_{k}}))\right]$
$\bm{\mathcal{C}}[c](\lambda_{k})\leftarrow$UniformRand($\mathcal{R}_{\lambda_{k}},1,1$)

43 end forall

44

45 end for

46 $\bm{\mathcal{C}}^{(*)}\leftarrow\operatorname*{arg\,min}\limits_{0\leq
j\leq b}\Psi(\bm{\mathcal{C}}[j])$ return
$\\{(\lambda_{i},\bm{\mathcal{C}}^{(*)},\Psi(\bm{\mathcal{C}}^{(*)})\\}$

47 end

To better understand the suggested collaborative randomized tuning process of
agents, an illustrative example is depicted in figure 2. In this figure, each
agent is represented by a different color, and the best candidate that each
agent finds at the end of each iteration is shown by a filled shape. Moreover,
we have assumed that the value of the loss function becomes smaller as we move
inwards in the depicted contour lines, and to prevent any exploration in the
domain of the subsidiary hyper-parameters, we have set
$\mathcal{E}=\\{\epsilon_{\lambda_{1}}=\frac{1}{6},\epsilon_{\lambda_{2}}=0\\}$
and
$\mathcal{E}=\\{\epsilon_{\lambda_{1}}=0,\epsilon_{\lambda_{2}}=\frac{1}{6}\\}$
for agents $g^{1}_{\lambda_{1}}$ and $g^{1}_{\lambda_{2}}$ respectively,
assuming that the domain size of each hyper-parameter is 1 and $b=3$. In
Iteration 1, both agents start at the top right corner of the search space and
are able to find candidates that yield lower loss function values than the
starting coordinate. For iteration 2, the starting coordinate of each agent is
set to the coordinate of the best candidate found by all agents in the
previous iteration. As the best candidate was found by agent
$g^{1}_{\lambda_{2}}$, we only see the change in the searching direction of
the red agent, i.e. $g^{1}_{\lambda_{2}}$. The winner agent at the end of this
iteration is agent $g^{1}_{\lambda_{1}}$, hence we do not see any change to
its searching direction in iteration 3. Please note that having 4 circles for
agent $g^{1}_{\lambda_{1}}$ in the last depicted iteration is because of
showing the starting coordinate, which happens to stay the best candidate, in
this iteration. It is also worth emphasizing that the starting coordinates are
not evaluated again by the agents, as they have already been accompanied by
their corresponding response values from the previous iterations.

Figure 2: A toy example demonstrating 3 iterations of running the proposed
method for tuning two hyper-parameters $\lambda_{1}$, and $\lambda_{2}$ using
terminal agents $g_{\lambda_{1}}^{1}$ and $g_{\lambda_{2}}^{1}$ respectively.
It is assumed that for each agent $b=3$.

##  3 Results and Discussion

This section dissects the performance of the proposed method in more detail.
It begins with the computational complexity of the technique, and then,
provides empirical results on both machine learning and general function
optimization tasks.

###  3.1 Computational Complexity

Forming the hierarchical structure and conducting the collaborative searching
process are the two major stages of the proposed method. These stages need to
be conducted in sequence, and the rest of this section investigates the
complexity of each step separately and in relation to one another.

Regarding the structural formation phase of the suggested method, the shape of
the hierarchy depends on the maximum number of connections that each agent can
handle; the fewer the number of manageable concurrent connections, the deeper
the resulting hierarchy. Using the same notations presented in section 2.1 and
assuming the same $c>1$ for all agents, the depth of the formed hierarchy is
$\lceil\log_{c}{|\bm{\lambda_{o}}|}\rceil$. Thanks to the distributed nature
of the formation algorithm and the concurrent execution of the agents, the
worst-case time complexity of the first stage will be
$\mathcal{O}(\log_{c}{|\bm{\lambda_{o}}|})$. With the same assumption, it can
be easily shown that the resulting hierarchical structure is a complete tree.
Hence, denoting the total number of agents in the system by $\mathbb{G}$, this
quantity would be:

$\frac{c^{\lceil\log_{c}{|\bm{\lambda_{o}}|}\rceil}-1}{c-1}<\mathbb{G}\leq\frac{c^{\lceil\log_{c}{|\bm{\lambda_{o}}|}\rceil+1}-1}{c-1}$
(6)

With that said, the space complexity for the first phase of the proposed
technique would be
$\mathcal{O}(\frac{c^{\lceil\log_{c}{|\bm{\lambda_{o}}|}\rceil+1}-1}{c-1})=\mathcal{O}(|\bm{\lambda_{o}}|)$.
It is worth noting that among all created agents, only $|\bm{\lambda_{o}}|$
terminal agents would require dedicated computational resources as they are
doing the actual searching and optimization process, and the remaining
$\mathbb{G}-|\bm{\lambda_{o}}|$ can all be hosted and managed together.

The procedures in each round of the second phase of the suggested method can
be broken into two main components: (i) transmitting the start coordinates
from the root of the hierarchy to the terminal agents, transmitting the
results back to the root, and preparing the feedback; and (ii) conducting the
actual searching process by the terminal agents to locate a local optimum. The
worst case time complexity of preparing the feedback based on the algorithms
that were discussed in section 2 would be $\mathcal{O}(|\bm{\lambda_{o}}|)$,
which is because of finding the best candidate among all returned results. In
addition, due to the concurrency of the agents, the first component is only
processed at the height of the built structure. Therefore, the time complexity
of component (i) would be
$\mathcal{O}(|\bm{\lambda_{o}}|+\log_{c}{|\bm{\lambda_{o}}|})=\mathcal{O}(|\bm{\lambda_{o}}|)$.
The complexity of the second component, on the other hand, depends on both the
budget of the agent, i.e. $b$, and the complexity of building and evaluating
the response function $\Psi$. Let $\mathcal{O}(\mathcal{R})$ denote the time
complexity of a single evaluation. As a terminal agent makes $b$ number of
such evaluations to choose its candidate optima, the time complexity for the
agent would be $\mathcal{O}(b\mathcal{R})$. As all agents work in parallel,
the complexity of a single iteration at the terminal agents would be
$\mathcal{O}(b\mathcal{R})$, leading to the overall time complexity of
$\mathcal{O}(|\bm{\lambda_{o}}|+b\mathcal{R})$. In machine learning problems,
we often have $\mathcal{O}(|\bm{\lambda_{o}}|)\ll\mathcal{O}(\mathcal{R})$.
Therefore, if $\mathcal{I}$ denotes the number of iterations until the second
phase of the tuning method stops, the complexity of the second stage would be
$\mathcal{O}(\mathcal{I}b\mathcal{R})$. The space complexity of the second
phase of the tuning method depends on the way that each agent is implementing
the main functionalities, such as the learning algorithms they represent,
transmitting the coordinates, and providing feedback. Except for the ML
algorithms, all internal functionalities of each agent can be implemented
using $\mathcal{O}(|\bm{\lambda_{o}}|)$ space. On the other hand, we have
$\mathbb{G}$ agents in the system, which leads to total space complexity of
$\mathcal{O}(|\bm{\lambda_{o}}|^{2})$ for non-ML tasks. Let
$\mathcal{O}(\mathcal{S})$ denote the worst-case space complexity of a machine
learning algorithm that we are tuning. The total space complexity of the
second phase of the proposed tuning method would be
$\mathcal{O}(|\bm{\lambda_{o}}|^{2}+\mathcal{S})$. Similar to the time
complexity, in machine learning, we often have
$\mathcal{O}(|\bm{\lambda_{o}}|)\ll\mathcal{O}(\mathcal{S})$, which makes the
total space complexity of second phase $\mathcal{O}(\mathcal{S})$. Please note
that we have factored out the budgets of the agents and the number of
iterations because of not storing the history between different evaluations
and iterations.

Considering both stages of the proposed technique, and due to the fact that
they are conducted in sequence, the time complexity of the entire steps in an
ML hyper-parameter tuning problem, from structure formation to completing the
searching operations, would be
$\mathcal{O}(\log_{c}{|\bm{\lambda_{o}}|}+\mathcal{I}b\mathcal{R})=\mathcal{O}(\mathcal{I}b\mathcal{R})$.
Similarly, the space complexity would be
$\mathcal{O}(|\bm{\lambda_{o}}|+\mathcal{S})=\mathcal{O}(\mathcal{S})$.

###  3.2 Empirical Results

This section presents the empirical results of employing the proposed agent-
based randomized searching algorithm and discusses the improvements resulting
from the suggested inter-agent collaborations. Hyper-parameter tuning in
machine learning is basically a back-box optimization problem, and hence, to
enrich our empirical discussions, this section also includes results from
multiple multi-dimensional optimization problems.

The performance metrics used for the experiments are based on the metrics that
are commonly used by the ML and optimization community. Additionally, we
analyze the behavior of the suggested methodology based on its own design
parameter values, such as budget, width, etc. The methods that have been
chosen for the sake of comparison are the standard random search, and the
Latin hypercube search methods [1] that are commonly used in practice. Our
choices are based on the fact that not only these methods are embarrassingly
parallel and are among the top choices to be considered in distributed
scenarios, but also they are used as the core optimization mechanisms of the
terminal agents in the suggested method and hence can present the impact of
the inter-agent collaborations better. In its generic format, as emphasized in
[21], one can easily employ alternative searching methods or diversify them at
the terminal level as needed.

Throughout the experiments, each terminal agent runs on a separate process,
and to make the comparisons fair we keep the number of model/function
evaluations fixed among all the experimented methods. To put it in more
detail, for a budget value of $b$ for each of $|\bm{\lambda_{o}}|$ terminal
agents and $\mathcal{I}$ number of iterations, the proposed method will
evaluate the search space in $b\times\mathcal{I}$ coordinates. We use the same
$|\bm{\lambda_{o}}|$ number of independent agents for the compared random-
based methodologies and keeping the evaluation budgets of the agents fixed —
the budgets are assumed to be enforced by the computational limitations of
devices or processes running the agents — we repeat those methods
$\mathcal{I}$ times and report the best performance among all agents’
repetition history as their final result.

The experiments assess the performance of the proposed method in comparison to
the other random-based techniques in 4 categories: (1) iteration-based
assessment, which checks the performance of the methods for a particular
iteration threshold. In this category, all other parameters, such as budget,
connection number, etc. are kept fixed; (2) budget-based assessment, which
examines the performance under various evaluation budgets for the terminal
agents. It is assumed that all agents have the same budget; (3) width-based
assessment, which checks how the proposed method performs for various
exploration criteria specified by the slot width parameter; and finally (4)
connection-based evaluation, which inspects the effect of the parallel
connection numbers that the internal agents can handle. In other words, this
evaluation checks if the proposed method is sensitive to the way that the
hyper-parameter or decision variables are split during the hierarchy formation
phase. All implementations use Python 3.9 and scikit-learn library [27], and
the results reported in all experiments are based on 50 different trials.

For the ML hyper-parameter tuning experiments, we have dissected the behavior
of the proposed algorithm in two classification and two regression problems.
The details of such problems, including the hyper-parameters that are tuned
and the used datasets, are presented in table 1. In all of the ML experiments,
we have used 5-fold cross-validation as the model evaluation method. The
results obtained for the classification and regression problems are plotted in
figures 3 and 4 respectively.

Table 1: The details of the machine learning algorithms and the datasets used
for hyper-parameter tuning experiments.

ML algorithm | $\bm{\lambda_{o}}$ | Dataset | Performance metric
---|---|---|---
C-Support Vector Classification (SVC) [28, 27] |  $\\{c,\gamma,\text{kernel}\\}$1 | artificial (100,20)$\dagger$ | accuracy
Stochastic Gradient Descent (SGD) Classifier [27] |  $\\{\alpha,\text{l1\\_ratio},\text{tol},\epsilon,\eta_{0},\text{val\\_frac}\\}$2 | artificial (500,20)$\dagger$ | accuracy
Passive Aggressive Regressor [29, 27] |  $\\{c,\text{tol},\epsilon,\text{val\\_frac}\\}$3 | artificial (300,100)$\ddagger$ | mean squared error
Elastic Net Regressor [30, 27] |  $\\{\alpha,\text{l1\\_ratio},\text{tol},\text{selection}\\}$4 | artificial (300,100)$\ddagger$ | mean squared error

  * 1

$c\sim logUniform(10^{-2},10^{13}),\gamma\sim
Uniform(0,1),\text{kernel}\in\\{\text{poly, linear, rbf, sigmoid}\\}$.

  * 2

$\alpha\sim Uniform(0,10^{3}),\text{l1\\_ratio}\sim
Uniform(0,1),\text{tolerance}\sim Uniform(0,10^{3}),\epsilon\sim
Uniform(0,10^{3}),\eta_{0}\sim Uniform(0,10^{3}),\\\
\text{validation\\_fraction}\sim Uniform(0,1)$.

  * 3

$c\sim Uniform(0,10^{3}),\text{tolerance}\sim
Uniform(0,10^{3}),\text{validation\\_fraction}\sim Uniform(0,1),\epsilon\sim
Uniform(0,1)$.

  * 4

$\alpha\sim Uniform(0,1),\text{l1\\_ratio}\sim
Uniform(0,1),\text{tolerance}\sim
Uniform(0,1),\text{selection}\in\\{cyclic,random\\}$.

  * $\dagger$

an artificially generated binary classification dataset using scikit-learn’s
make_classification function [31]. The first number represents the number of
samples and the second figure is the number of features.

  * $\ddagger$

an artificially generated regression dataset using scikit-learn’s
make_regression function [31]. The first number represents the number of
samples and the second figure is the number of features.

Figure 3: Average performance of C-Support Vector Classification (SVC) (first
row) and Stochastic Gradient Descent (SGD) (second row) classifiers on two
synthetic classification datasets based on the accuracy measure. The error
bars in each plot are calculated based on the standard error. Figure 4:
Average performance of Passive Aggressive (first row) and Elastic Net (second
row) regression algorithms on two synthetic regression datasets based on the
mean squared error (MSE) measure. The error bars in each plot are calculated
based on the standard error.

For the _iterations_ plot in the first column plots of figures 3 and 4, we
fixed the parameters of the proposed method for all agents as follows:
$b=3,\mathcal{E}=2^{-6},c=2,\Delta=\\{2,2,\dots,2\\}$. As can be seen, when
the proposed method is allowed to run for more iterations, it yields better
performance and its superiority against the other two random-based methods is
evident. Comparing the relative performance improvements resulting from the
proposed method in the presented ML tasks, it can be seen that as the search
space of the agents, and the number of the hyper-parameters needed to be tuned
increase, the proposed collaborative method achieved higher improvement. For
the SGD classifier, for instance, the objective hyper-parameter set comprises
6 members with continuous domain spaces, and the amount of improvements that
have been made after 10 iterations is much higher, about 17%, than the other
experiments with 3-4 hyper-parameters and mixed continuous and discrete domain
spaces.

The second column of figures 3 and 4 illustrates how the performance of the
proposed technique changes when we increase the evaluation budgets of the
terminal agents. For this set of experiments, we set the parameter values of
our method as follows:
$\mathcal{I}=10,\mathcal{E}=2^{-6},c=2,\Delta=\\{2,2,\dots,2\\}$. By
increasing the budget value, the performance of the suggested approach per se
improves. However, the rate of improvement slows down for higher budget
values, and comparing it against the performance of the other two random-based
searching methods, the improvement is significant for lower values of budget.
In other words, the proposed tuning method surpasses the other two methods
when the agents have limited searching resources. This makes our method a good
candidate for tuning the hyper-parameters of deep learning approaches with
expensive model evaluations.

The behavior of the suggested method under various exploration parameter
values can be seen in column 3 of figures 3 and 4. The $\omega$ values on the
x-axis of the plots are used to set the initial value for the slot width
parameter of all agents using $\mathcal{E}=2^{-\omega-1}$. Based on this
configuration, higher values of $\omega$ yield lower values of $\mathcal{E}$,
and as a result, more exploitation around the starting coordinates. The other
parameters of the method are configured as follows:
$\mathcal{I}=10,b=3,c=2,\Delta=\\{2,2,\dots,2\\}$. Recall from section 2 that
the exploration parameter is used by an agent for the dimensions that it does
not represent. Based on the results obtained from various tasks, choosing a
proper value for this parameter depends on the characteristics of the response
function. Having said that, the behavior for a particular task remains almost
consistent. Hence, trying one small and one large value for this parameter in
a specific problem will reveal its sensitivity and help choose an appropriate
value for it.

Finally, the last set of experiments investigates the impact of the number of
parallel connections that the internal agents can manage, i.e. $c$, on the
performance of the suggested method. The results of this study are plotted on
the last column of figures 3 and 4. The difference in the number of data
points in each plot is because of the difference in the size of the hyper-
parameters that we tune in each task. The values of the parameters that we
kept fixed for this set of experiments are as follows:
$\mathcal{I}=10,b=3,\mathcal{E}=2^{-6},\Delta=\\{2,2,\dots,2\\}$. As can be
seen from the illustrated results, the proposed method is not very sensitive
to the value that we choose or is enforced by the system for parameter $c$.
This parameter plays a critical role in the shape of the hierarchy that is
distributedly formed in phase 1 of the suggested approach; therefore, one can
opt to choose a value that fits with the connection or computational resources
that are available without sacrificing the performance very much.

As stated before, we have also studied the suggested technique for the global
optimization problem to see how it performs in finding the optima of various
convex and non-convex functions. These experiments also help us closely check
the relative performance improvements in higher dimensions. We have chosen
three non-convex benchmark optimization functions and a convex toy function
the details of which are presented in table 2. For each function, we run the
experiments in three different dimension sizes, and the goal of the
optimization is to find the global minimum. Very similar to the settings that
we discussed for ML hyper-parameter tuning, whenever we mean to fix the value
of each parameter values in different experiment sets, we use the following
parameter values:
$\mathcal{I}=10,b=3,\mathcal{E}=2^{-10},c=2,\Delta=\\{2,2,\dots,2\\}$.

Table 2: The details of the multi-dimensional functions used for global
optimization experiments.

Function | $\bm{\lambda_{o}}$ | Domain | $f(\bm{x}^{*})$
---|---|---|---
Hartmann, 3D, 4D, 6D [32] | $\\{x_{1},\dots,x_{d}\\},d\in\\{3,4,6\\}$ | $x_{i}\in\left[0,1\right]$ |  3D:-3.86278, 4D:-3.135474, 6D:-3.32237
Rastrigin, 3D, 6D, 10D [33] | $\\{x_{1},\dots,x_{d}\\},d\in\\{3,6,10\\}$ | $x_{i}\in\left[-5.12,5.12\right]$ |  3D:0,6D:0, 10D:0
Styblinski-Tang, 3D, 6D, 10D [34] | $\\{x_{1},\dots,x_{d}\\},d\in\\{3,6,10\\}$ | $x_{i}\in\left[-5,5\right]$ |  3D:-117.4979,6D:-234.9959, 10D:391.6599
Mean Average Error, 3D, 6D, 10D$\dagger$ | $\\{x_{1},\dots,x_{d}\\},d\in\\{3,6,10\\}$ | $x_{i}\in\left[0,100\right]$ |  3D:0,6D:0, 10D:0$\ddagger$

  * $\dagger$

This is a toy multi-dimensional MAE function that is defined as
$f(\bm{x})=\frac{1}{n}\sum_{i=1}^{n}{|\bm{x}-\chi|}$ where $\chi$ denotes a
ground truth vector that is generated randomly in the domain space for each
experiment.

  * $\ddagger$

This is a convex function and the coordinate of its minimum value depends on
the ground truth vector that is generated, i.e. when $\bm{x}=\chi$.

The plots are grouped by functions and can be found in figures 5, 6, 7, and 8.
The conclusion that was made about the behavior of the proposed approach under
different values of its design parameters applies in these optimization
experiments as well. That is, the more the proposed method runs, the better
performance it achieves; its superiority on low budget values is clear; its
sensitivity to exploration parameter values is consistent; and the way that
the decision variables set is broken down during the formation of the
hierarchy does not affect the performance very much. Furthermore, as can be
seen in each group figure, the proposed algorithm yields a better minimum
point in comparison to the other two random-based methods when the
dimensionality of a function increases.

Figure 5: Average function values of the Hartmann function under variable
iteration, budget, exploration, and connection thresholds. Each row of the
figure pertains to a particular dimension size and the error bars are
calculated based on the standard error. Figure 6: Average function values of
the Rastrigin function under variable iteration, budget, exploration, and
connection thresholds. Each row of the figure pertains to a particular
dimension size and the error bars are calculated based on the standard error.
Figure 7: Average function values of the Styblinski-Tang function under
variable iteration, budget, exploration, and connection thresholds. Each row
of the figure pertains to a particular dimension size and the error bars are
calculated based on the standard error. Figure 8: Average function values of
the Toy Mean Absolute Error function under variable iteration, budget,
exploration, and connection thresholds. Each row of the figure pertains to a
particular dimension size and the error bars are calculated based on the
standard error.

It is worth emphasizing that the contribution of this paper is not to compete
with the state-of-the-art in function optimization but to propose a
distributed tuning/optimization approach that can be deployed on a set of
distributed and networked devices. The discussed analytical and empirical
results not only demonstrated the behavior and impact of the design parameters
that we have used in our approach but also suggested the way that they can be
adjusted for different needs. We believe the contributed approach of this
paper can be significantly improved provided with more sophisticated and
carefully chosen tuning strategies and corresponding configurations.

##  4 Conclusion

This paper presented an agent-based collaborative random search method that
can be used for machine learning hyper-parameter tuning and global
optimization problems. The approach employs two types of agents during the
tuning/optimization process: the internals and the terminal agents that are
responsible for facilitating the collaborations and tuning an individual
decision variable respectively. Such agents and the interaction network
between them are created during the hierarchy formation phase and remain the
same for the entire run-time of the suggested method. Thanks to the modular
and distributed nature of the approach and its procedures, it can be easily
deployed on a network of devices with various computational capabilities.
Furthermore, the design parameters used in this technique enable each
individual agent to customize its own searching process and behavior
independent from its peers in the hierarchy, and hence, allow employing
diversity in both algorithmic and deployment levels.

The paper dissected the proposed model from different aspects and provided
some tips on handling its behavior for various applications. According to the
analytical discussions, our approach requires slightly more amount of
computational and storage resources than the traditional and Latin hyper-cube
randomized search methods that are commonly used for both hyper-parameter
tuning and black-box optimization problems. However, this results in
significant performance improvements, especially in computationally restricted
circumstances and higher numbers of decision variables. This conclusion was
verified in both machine learning model tuning tasks and general multi-
dimensional function optimization problems.

The presented work can be further extended both technically and empirically.
As was discussed throughout this paper, we kept the searching strategies and
the way the design parameters are configured as simple as possible so we could
reach a better understanding of the effectiveness of collaborations and
searching space divisions. A few potential extensions in this direction
include: utilizing diverse searching methods and hence having a heterogeneous
multi-agent system at the terminal level; splitting the searching space not
based on the dimensions but based on the range of the values that decision
variables in each dimension can have; employing more sophisticated
collaboration techniques; and using a learning-based approach to dynamically
adapt the values of the design parameters during the runtime of the method.
Empirically, on the other hand, the presented research can be extended by
applying it to expensive machine learning tasks such as tuning deep learning
models with a large number of hyper-parameters. We are currently working on
some of these studies and suggest them as future work.

## References

  * [1] James Bergstra and Yoshua Bengio.  Random search for hyper-parameter optimization.  Journal of machine learning research, 13(2), 2012.
  * [2] Ron Kohavi and George H John.  Automatic parameter selection by minimizing estimated error.  In Machine Learning Proceedings 1995, pages 304–312. Elsevier, 1995\.
  * [3] B. Bischl, O. Mersmann, H. Trautmann, and C. Weihs.  Resampling Methods for Meta-Model Validation with Recommendations for Evolutionary Computation.  Evolutionary Computation, 20(2):249–275, 06 2012.
  * [4] Douglas C Montgomery.  Design and analysis of experiments.  John wiley & sons, 2017.
  * [5] George H. John.  Cross-validated c4.5: Using error estimation for automatic parameter selection.  Technical report, 1994.
  * [6] Jonas Močkus.  On bayesian methods for seeking the extremum.  In Optimization techniques IFIP technical conference, pages 400–404. Springer, 1975.
  * [7] Jonas Mockus.  Bayesian approach to global optimization: theory and applications, volume 37.  Springer Science & Business Media, 2012.
  * [8] Matthias Feurer and Frank Hutter.  Hyperparameter optimization.  In Automated machine learning, pages 3–33. Springer, Cham, 2019\.
  * [9] Dan Simon.  Evolutionary optimization algorithms.  John Wiley & Sons, 2013.
  * [10] A Esmaeili and N Mozayani.  Adjusting the parameters of radial basis function networks using particle swarm optimization.  In 2009 IEEE International Conference on Computational Intelligence for Measurement Systems and Applications, pages 179–181. IEEE, 2009\.
  * [11] Hussain Alibrahim and Simone A. Ludwig.  Hyperparameter optimization: Comparing genetic algorithm against grid search and bayesian optimization.  In 2021 IEEE Congress on Evolutionary Computation (CEC), pages 1551–1559, 2021.
  * [12] Richard E Bellman.  Adaptive control processes.  Princeton university press, 1961.
  * [13] Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P Adams, and Nando De Freitas.  Taking the human out of the loop: A review of bayesian optimization.  Proceedings of the IEEE, 104(1):148–175, 2015.
  * [14] Javier Garcia-Barcos and Ruben Martinez-Cantin.  Fully distributed bayesian optimization with stochastic policies.  2019\.
  * [15] M Todd Young, Jacob D Hinkle, Ramakrishnan Kannan, and Arvind Ramanathan.  Distributed bayesian optimization of deep reinforcement learning algorithms.  Journal of Parallel and Distributed Computing, 139:43–52, 2020\.
  * [16] Peter I. Frazier.  A tutorial on bayesian optimization, 2018.
  * [17] Frauke Friedrichs and Christian Igel.  Evolutionary tuning of multiple svm parameters.  Neurocomputing, 64:107–117, 2005.
  * [18] Ilya Loshchilov and Frank Hutter.  Cma-es for hyperparameter optimization of deep neural networks.  arXiv preprint arXiv:1604.07269, 2016.
  * [19] Dominik Ryzko.  Modern big data architectures: a multi-agent systems perspective.  John Wiley & Sons, 2020.
  * [20] Ahmad Esmaeili, John C. Gallagher, John A. Springer, and Eric T. Matson.  Hamlet: A hierarchical agent-based machine learning platform.  ACM Transactions on Autonomous and Adaptive Systems, 16(3–4), 2022\.
  * [21] Ahmad Esmaeili, Zahra Ghorrati, and Eric T. Matson.  Hierarchical collaborative hyper-parameter tuning.  In Advances in Practical Applications of Agents, Multi-Agent Systems, and Complex Systems Simulation. The PAAMS Collection, pages 127–139, Cham, 2022. Springer International Publishing.
  * [22] Rémi Bardenet, Mátyás Brendel, Balázs Kégl, and Michele Sebag.  Collaborative hyperparameter tuning.  In International conference on machine learning, pages 199–207. PMLR, 2013.
  * [23] Thomas Swearingen, Will Drevo, Bennett Cyphers, Alfredo Cuesta-Infante, Arun Ross, and Kalyan Veeramachaneni.  Atm: A distributed, collaborative, scalable system for automated machine learning.  In 2017 IEEE international conference on big data (big data), pages 151–162. IEEE, 2017.
  * [24] Patrick Koch, Oleg Golovidov, Steven Gardner, Brett Wujek, Joshua Griffin, and Yan Xu.  Autotune: A derivative-free optimization framework for hyperparameter tuning.  In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 443–452, 2018.
  * [25] Arman Iranfar, Marina Zapater, and David Atienza.  Multi-agent reinforcement learning for hyperparameter optimization of convolutional neural networks.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 2021.
  * [26] Jack Parker-Holder, Vu Nguyen, and Stephen J Roberts.  Provably efficient online hyperparameter optimization with population-based bandits.  Advances in Neural Information Processing Systems, 33:17200–17211, 2020.
  * [27] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [28] Chih-Chung Chang and Chih-Jen Lin.  Libsvm: a library for support vector machines.  ACM transactions on intelligent systems and technology (TIST), 2(3):1–27, 2011.
  * [29] Koby Crammer, Ofer Dekel, Joseph Keshet, Shai Shalev-Shwartz, and Yoram Singer.  Online passive aggressive algorithms.  Journal of Machine Learning Research, 7:551–585, 2006.
  * [30] Jerome Friedman, Trevor Hastie, and Rob Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  Journal of statistical software, 33(1):1, 2010.
  * [31] Scikit-learn api reference.  https://scikit-learn.org/stable/modules/classes.html.  Accessed: 2022-11-24.
  * [32] Momin Jamil and Xin She Yang.  A literature survey of benchmark functions for global optimisation problems.  International Journal of Mathematical Modelling and Numerical Optimisation, 4(2):150, 2013.
  * [33] Günter Rudolph.  Globale Optimierung mit parallelen Evolutionsstrategien.  PhD thesis, Diplomarbeit, Universit at Dortmund, Fachbereich Informatik, 1990.
  * [34] MA Styblinski and T-S Tang.  Experiments in nonconvex optimization: stochastic approximation with function smoothing and simulated annealing.  Neural Networks, 3(4):467–483, 1990.


# Evidence for strong electron correlations in a non-symmorphic Dirac
semimetal

Yu-Te Hsu<EMAIL_ADDRESS>High Field Magnet Laboratory (HFML-EMFL) and
Institute for Molecules and Materials, Radboud University, Toernooiveld 7,
6525 ED Nijmegen, Netherlands  Danil Prishchenko  Department of Theoretical
Physics and Applied Mathematics, Ural Federal University, 620002 Ekaterinburg,
Russia  Maarten Berben  Matija Čulo  Steffen Wiedmann  High Field Magnet
Laboratory (HFML-EMFL) and Institute for Molecules and Materials, Radboud
University, Toernooiveld 7, 6525 ED Nijmegen, Netherlands  Emily C. Hunter
School of Physics and Astronomy, The University of Edinburgh, James Clerk
Maxwell Building, Mayfield Road, Edinburgh EH9 2TT, United Kingdom  Paul
Tinnemans  Department of Solid State Chemistry, Radboud University,
Heyendaalseweg 135, 6525 AJ Nijmegen  Tomohiro Takayama  Max Planck Institute
for Solid State Research, Heisenbergstrasse 1, 70569 Stuttgart, Germany
Vladimir Mazurenko  Department of Theoretical Physics and Applied Mathematics,
Ural Federal University, 620002 Ekaterinburg, Russia  Nigel E. Hussey
<EMAIL_ADDRESS>High Field Magnet Laboratory (HFML-EMFL) and Institute for
Molecules and Materials, Radboud University, Toernooiveld 7, 6525 ED Nijmegen,
Netherlands  H. H. Wills Physics Laboratory, University of Bristol, Tyndall
Avenue, Bristol BS8 1TL, United Kingdom  Robin S. Perry<EMAIL_ADDRESS>London Centre for Nanotechnology and Department of Physics and Astronomy,
University College London, London WC1E 6BT, United Kingdom  ISIS Neutron and
Muon Source, Rutherford Appleton Laboratory, Harwell, OX11 0QX, United Kingdom

###### Abstract

Metallic iridium oxides (iridates) provide a fertile playground to explore new
phenomena resulting from the interplay between topological protection, spin-
orbit and electron-electron interactions. To date, however, few studies of the
low energy electronic excitations exist due to the difficulty in synthesising
crystals with sufficiently large carrier mean-free-paths. Here, we report the
observation of Shubnikov-de Haas quantum oscillations in high-quality single
crystals of monoclinic SrIrO3 in magnetic fields up to 35 T. Analysis of the
oscillations reveals a Fermi surface comprising multiple small pockets with
effective masses up to 4.5 times larger than the calculated band mass. Ab-
initio calculations reveal robust linear band-crossings at the Brillouin zone
boundary, due to its non-symmorphic symmetry, and overall we find good
agreement between the angular dependence of the oscillations and the
theoretical expectations. Further evidence of strong electron correlations is
realized through the observation of signatures of non-Fermi liquid transport
as well as a large Kadowaki-Woods ratio. These collective findings, coupled
with knowledge of the evolution of the electronic state across the Ruddlesden-
Popper iridate series, establishes monoclinic SrIrO3 as a topological
semimetal on the boundary of the Mott metal-insulator transition.

## Introduction

Band topology and strong electron correlations represent two of the most
active research themes in quantum materials [1], with much attention now
focused on the search for physics arising from their coexistence. To date,
however, few examples of correlated topological materials are known to exist
[2, 3, 4, 5]. Among these, iridium oxides have emerged as one of the most
promising material platforms on which to investigate the interplay between
spin-orbit and electron-electron interactions [6, 7, 8], due to the comparable
energy scales of spin-orbit coupling ($\Lambda$), Coulomb repulsion ($U$) and
electron bandwidth ($W$) [9].

The Ruddlesden-Popper series Srn+1IrnO3n+1, for example, exhibits a plethora
of intriguing properties, such as spin-orbit-coupled Mottness [10], pseudogap
phenomenology [11, 12, 13, 14], odd-parity hidden order [15], and metal-
insulator transitions [16, 17, 18, 19]. The spin-orbit-coupled $J_{\rm
eff}=1/2$ and $J_{\rm eff}=3/2$ bands, resulting from crystal-field splitting
of the $e_{\rm g}$ and $t_{\rm 2g}$ orbitals, are often treated as the
starting point for their understanding [7]. In single-layer Sr2IrO4 ($n=1$),
the electron bandwidth of the half-filled $J_{\rm eff}=1/2$ band is such that
$U/W>1$, resulting in a spin-orbit-coupled Mott insulator. With increasing
$n$, electron hopping (and thereby $W$) increases while at the same time $U$
is reduced. Consequently, bilayer Sr3Ir2O7 is found to be only weakly
insulating [16] while the infinite-layer end member SrIrO3 is itinerant,
though its semimetallic ground state and low carrier concentration appear at
odds with theoretical expectations of a metallic state with a half-filled
$J_{\rm eff}=1/2$ band and a large Fermi surface. Recently, the
semimetallicity of SrIrO3 was proposed to have a topological origin [20, 21,
22].

SrIrO3 crystallizes in two polymorphs: an ambient-pressure monoclinic phase
and a high-pressure orthorhombic phase, which is also stabilized under
epitaxial strain [17, 18]. Three-dimensional (3D) Dirac points are created at
the zone boundary of the non-symmorphic lattice by glide symmetry and are
protected against gapping caused by spin-orbit coupling. This feature, coupled
with their perceived proximity to a Mott transition, makes them attractive
candidates for realizing a correlated topological semimetallic state. Until
now, however, studies of their intrinsic electronic structure and properties
have been restricted due to the fact that orthorhombic SrIrO3 ($o$-SIO3) can
only be synthesized in thin-film form at ambient pressure [9, 17, 18, 23, 19,
24, 25]), whereas monoclinic SrIrO3 ($m$-SIO3), which can be synthesized in
bulk form, has thus far remained largely unexplored.

Here, we report the determination of the full Fermi surface of $m$-SIO3 via
the observation of Shubnikov-de Haas oscillations in high-quality single
crystals. A number of small pockets are detected, in excellent agreement with
first-principles density-functional theory (DFT) calculations. According to
DFT, the electronic band structure of $m$-SIO3 is found to be highly sensitive
to atomic displacements, though the presence of linear (Dirac) band-crossings
at the A- and M-points in the first Brillouin zone remains robust.
Significantly, the experimentally determined effective masses are found to be
substantially enhanced compared to the DFT band mass, signifying the presence
of strong electron correlations. We find a good overall agreement between the
theoretically derived angular dependence of the quantum oscillation
frequencies with the experimental observation, indicating that correlation
effects play only a minor role in determining the Fermi surface topology of
$m$-SIO3. Further evidence for strong correlation effects is provided by the
observation of a linear-in-temperature ($T$) component in the low-$T$
resistivity as well as an unusually large Kadowaki-Woods ratio. Overall, this
study demonstrates that $m$-SIO3 is a rare iridate topological semimetal
displaying signatures of both a correlated Fermi liquid and a non-Fermi liquid
ground state, possibly linked to its proximity to a Mott transition. The
availability of high-quality single crystals and the sensitivity of its
electronic structure to small perturbations make $m$-SIO3 a promising platform
on which to explore correlation-driven physics in a topological system.

## Results and Discussions

### Shubnikov-de Haas oscillations

Figure 1a shows field-dependent resistivity $\rho_{xx}(B)$ curves measured at
$T$ = 0.36 K and two angles $\theta=70^{\circ}$ and 83∘ relative to the
[001]-axis (see inset to Fig. 1a for experimental alignment). Prominent
Shubnikov-de Haas (SdH) oscillations emerge above 20 T (most evident in the
$\theta=83^{\circ}$ trace). By plotting the field-derivative d$\rho_{xx}$/d$B$
in Fig. 1b, the full oscillatory signals become visible for both angles. The
oscillation waveform at $\theta=83^{\circ}$, which consists of four
oscillatory components as revealed by the fast Fourier transform (FFT)
spectrum, is dominated by the $F_{\rm SdH}$ = 650 T component. At
$\theta=70^{\circ}$, all four constituent frequencies, with distinct peaks at
$F_{\rm SdH}$ = 75, 422, 653, and 907 T, can be tracked over a range of
temperatures, as shown in Fig. 1c, thereby enabling a more reliable
determination of the corresponding quasiparticle effective masses $m^{*}$.
Each frequency $F_{\rm SdH}$ is related to an extremal cross-sectional area
($A_{k}$) of the Fermi surface perpendicular to the field direction via the
Onsager relation $F_{\rm SdH}=(\hbar/2\rm{\pi e})$$A_{k}$. By fitting the
$T$-dependence of the oscillation amplitude $A_{\rm FFT}(T)$ at
$\theta=70^{\circ}$ using the Lifshitz-Kosevich expression for the thermal
damping term [26] $R_{T}=X/\sinh X,$ where $X$ = 14.69 $m^{*}T/B$, we obtain
$m^{*}$ values associated with each frequency ranging from 1.9 to 4.8 $m_{\rm
e}$. Such high masses are exceptional for a Dirac semimetal.

Figure 1: Shubnikov-de Haas oscillations in monoclinic SrIrO3. a Electrical
resistivity $\rho_{xx}$ measured in a magnetic field B up to 35 T at specified
angles. Inset: crystal structure (in conventional unit cell) and experimental
alignment. B is rotated within the $a$-$c$ plane, with $\theta$ denoting the
angle between B and the $c$-axis [001]. b Field-derivative of the data shown
in a using the same colour coding. $\theta=70^{\circ}$ data are shifted
vertically for clarity. c Fourier spectrum of the SdH oscillations measured at
$\theta=70^{\circ}$ and indicated temperatures, after a polynomial background
is subtracted. Four distinct peaks at $F_{\rm SdH}$ = 75, 422, 653, and 907 T
are resolved using a field window of 18 T $<B<$ 35 T. d Extraction of the
effective masses $m^{*}$ through fitting the oscillation amplitude $A_{\rm
FFT}(T)$ using the Lifshitz-Kosevich expression (see text). Error bars in
$A_{\rm FFT}$ reflect the noise floor of the Fourier spectra shown in c.

### DFT calculations

Figure 2: Electronic structure of $m$-SIO3 calculated by density-functional
theory. a, b Band structures obtained using room-temperature (RT) lattice
parameters reported in [27] and low-temperature (LT) lattice parameters found
at 13 K on our crystals. The high-symmetry points in the first Brillouin zone
are labelled in f. c, d Comparison of the RT and LT band structures with and
without SOC near the Fermi energy $\epsilon_{\rm F}$. Bands that give rise to
electron and hole Fermi pockets are highlighted in cyan and yellow,
respectively. Note that linear band-crossings at the M- and A-points are found
in both GGA+SOC results despite the overall differences. e, f Corresponding
Fermi surface models for the RT and LT band structures (GGA+SOC). Fermi
pockets are colour-coded in accordance with c, d.

In order to understand the origin of the observed SdH oscillations, we
performed DFT band structure calculations using the generalized gradient
approximation with effects of spin-orbit coupling included (GGA+SOC). Figure 2
shows the calculated band structures using the structural parameters reported
in [27], refined at room temperature (RT), and those found on our own single
crystals, refined at low temperature (LT; 13 K) (see Supplementary Methods).
The RT band structure agrees well with a previous report [22], with linear
band-crossings near the Fermi energy $\epsilon_{\rm F}$ at the M- and
A-points. The LT band structure, however, reveals a markedly different
picture, though notably the linear band-crossings at M- and A-points are still
preserved.

The effects of SOC on each band structure are illustrated in Fig. 2c,d. The
inclusion of SOC leads to a separation of the low-lying bands near
$\epsilon_{\rm F}$ by up to $\sim$ 0.3 eV, while the degeneracy at the M- and
A-points remains unaffected. The resultant Fermi surfaces from the two GGA+SOC
band structures, shown in Fig. 2e,f, show a number of notable differences
which we attribute to the extreme sensitivity of the electronic structure of
$m$-SIO3 with respect to the precise atomic positions – known to be a
challenging issue for iridates [25, 28] – rather than to a structural phase
transition. The different types of samples used for the structural refinements
– polycrystalline powder in [27] and crushed single crystals in our study –
can cause minute structural differences and, in turn, lead to differences in
the corresponding band structures. Nevertheless, the linear band-crossings
near $\epsilon_{\rm F}$ at the M- and A-points are found to be a robust
feature independent of such structural details, confirming their protection by
the underlying non-symmorphic lattice symmetry.

Figure 3: Angular dependence of the observed and calculated quantum
oscillation frequencies. a $F_{\rm SdH}(\theta)$ measured in two samples, S1
and S2, cut from different crystals of the same growth batch. Four distinct
branches can be identified, labelled as $F_{1}$ to $F_{4}$. b, c $F_{\rm
SdH}(\theta)$ expected from the LT and RT Fermi surface models as shown in
Fig. 2, with the electron and hole-pockets colour-coded accordingly. Note that
the $y$-axis for all panels is the same.

In order to identify the appropriate Fermi surface model from the DFT
calculations, we compare in Fig. 3 the angular dependence of the observed
quantum oscillation frequencies $F_{\rm SdH}(\theta)$ with the predictions
emerging from each model. As shown in Fig. 3a, the four distinct frequency
branches disperse differently with angle $\theta$ and range from 50 T to 2000
T. Overall, we find a better agreement between the experiment with the LT-SOC
model (referred to as the LT model hereafter), which reproduces the close
tracking between $F_{1}$ and $F_{2}$, the rapid increase in $F_{\rm SdH}$ as
$\theta\rightarrow 0$, as well as the moderate variation of $F_{3}$ with
respect to $\theta$. This observation implies all the observed pockets
originate from the metallic bulk. While the apparent divergence of frequencies
$F_{1,2}$ as B approaches [001] suggests a possible surface state originating
in the (100) plane, we can exclude such a possibility by examining its precise
angle dependence. For a surface orbit, the SdH frequency is expected to follow
a $1/\cos\phi$ dependence as B is tilted at an angle $\phi$ away from the
surface normal in which the orbit is hosted. As shown in Supplementary Fig. 3,
by plotting $F_{1,2}\cos\phi$ as a function of $\phi$, where
$\phi=90^{\circ}-\theta$ denotes the tilt angle between B and [100], we find
$F_{1,2}\cos\phi$ deviates strongly from a constant value as $\phi$ exceeds
30∘, inconsistent with expectations for a surface origin. We therefore
conclude that there is no surface contribution to the observed SdH
frequencies. We note that the LT model underestimates the $F_{\rm SdH}$ for
the $F_{1,2}$ branch as $\theta\rightarrow$ 90∘, indicating the electron
pocket at M-point is less elongated in reality than that illustrated in Fig.
2f, possibly due to subtle atomic displacements and/or correlation effects.

Nevertheless, despite the technical challenges associated with DFT
calculations for this complex oxide, the overall agreement between the
experimentally observed SdH frequencies and the theoretical expectations from
the LT model suggests that the full Fermi surface of $m$-SIO3 has been
determined by our study. The Dirac points lie $\approx$ 15 and 50 meV below
the Fermi level at the M- and A-points, respectively, at which the band
dispersion remains linear (see Fig. 2d). In order to identify the correct FS
model from the DFT calculations, we compare in Fig. 3 the angular dependence
of the observed SdH frequencies $F_{\rm SdH}(\theta)$ with the predictions
emerging from each model. As shown in Fig. 3(a), the four distinct frequency
branches disperse differently with angle $\theta$ and range from 50 T to 2000
T. While both models predict four branches at most angles and an electron-
pocket associated with the highest frequency branch (Fig. 3(b, c)), we find
better overall agreement with the LT-SOC model (referred to as the LT model
hereafter). Importantly, the close tracking between $F_{1}$ and $F_{2}$, the
rapid increase in $F_{\rm SdH}$ as $\theta\rightarrow 0$, and the moderate
variation of $F_{3}$ with respect to $\theta$ are all well-reproduced in the
LT model.

The rapid increase of $F_{1,2}$ as B is tilted away from [100] may suggest a
surface state located within the (100) plane. For such an orbit, the SdH
frequency should follow a $1/\cos\phi$ dependence as B is tilted away from the
surface. We examine this possibility by plotting $F_{1,2}\cos\phi$ versus
$\phi$, where $\phi=90^{\circ}-\theta$ denotes the tilt angle between B and
[100], as shown in Supplemental Material Fig. S2. We find $F_{1,2}\cos\phi$
deviates strongly from a constant value as $\phi$ exceeds 30∘, inconsistent
with expectations for a surface state. We therefore conclude that all the
observed pockets arise from the metallic bulk. Despite the technical
challenges associated with DFT calculations for this complex oxide, the
remarkable agreement in the experimentally observed angular dependence of SdH
frequencies and the theoretical expectations from the LT model suggests that
the full FS of $m$-SIO3 has been accurately determined by our study.

### Signatures of strong electronic correlations

Table 1: Fermi surface parameters found by experiment and DFT calculations (LT model). $l$, $m^{*}$, and $m_{\rm b}$ represent quasiparticle mean free path, effective mass, and DFT band mass, respectively. All parameters are extracted at $\theta=70^{\circ}$. Note that a reliable extraction of $l$ associated with $F_{4}$ is precluded by the dominating influence of the higher-frequency oscillations, and the parameters marked by $\dagger$ are an average over the multiple small orbits ($<$ 30 T) predicted by the DFT model that cannot be resolved experimentally. |  $F_{\rm SdH}$ (T) |  $F_{\rm DFT}$ (T) |  $l$ (nm) |  $m^{*}$ ($m_{\rm e}$) | $m_{\rm b}(m_{\rm e})$ | $m^{*}/m_{\rm b}$
---|---|---|---|---|---|---
$F_{1}$ | 907 | 462 | 30 $\pm$ 2 | 2.81 $\pm$ 0.03 | 1.46 | 1.9
$F_{2}$ | 653 | 410 | 39 $\pm$ 3 | 1.90 $\pm$ 0.03 | 1.36 | 1.4
$F_{3}$ | 422 | 318 | 26 $\pm$ 3 | 5.8 $\pm$ 0.7 | 1.07 | 5.4
$F_{4}$ | 75 | 15† | - | 2.0 $\pm$ 0.4 | 1.04† | 1.9†
Figure 4: Non-Fermi liquid transport and Kadowaki-Woods ratio of $m$-SIO3. a
Zero-field resistivity $\rho_{xx}(T)$ below 15 K and, b its derivative
d$\rho_{xx}$/d$T$ below 6 K, revealing a finite $T$-linear component in
$\rho_{xx}(T)$ at low $T$. Purple dashed line in a is a fit to the
experimental data using the functional form:
$\rho_{xx}(T)=\rho_{0}+A_{1}T+A_{2}T^{2}$, with the $T$-linear and
$T^{2}$-components shown in red and blue, respectively. c Specific heat
plotted as $C/T$ versus $T^{2}$ and fitted with $C/T=\gamma_{\rm 0}+\beta
T^{2}$ below 10 K (blue line). d Magnetic susceptibility $\chi_{\rm m}$
measured with $B$ = 2 T applied along [001]. $\chi_{\rm m}$ is largely
$T$-independent down to $\approx$ 30K, with a constant $\chi_{0}$ as marked by
the grey band. Below 30 K, $\chi_{\rm m}$ shows a small upturn, which we
attribute to an impurity contribution that follows a Curie-Weiss behaviour. e
Transport coefficient $A_{2}$ plotted against $\gamma_{0}^{2}$, known as the
Kadowaki-Woods plot, for selected correlated oxides (diamonds), iron pnictides
(hexagons), heavy fermion materials (circles) ([30, 31, 32] and references
therein), and monoclinic SrIrO3 (elongated diamond). Red line corresponds to
$A_{2}/\gamma_{0}^{2}$ = 10 $\mu\Omega$ cm mol2 K2 mJ-6, known to describe
many heavy fermion materials well.

Having established the appropriate Fermi surface model, we proceed to quantify
the extent of electron correlation in $m$-SIO3 by comparing the measured
effective masses $m^{*}$ with the calculated band masses $m_{\rm b}$. Here we
use the definition $m_{\rm b}=\hbar^{2}(dA_{k}/dE)|_{\epsilon_{\rm F}}$,
namely the change of area enclosed by the cyclotron orbit in $k$-space with
respect to the change in energy at the Fermi level, to describe the band mass
of quasiparticles with an arbitrary dispersion $\epsilon(k)$. As shown in
Table 1, we find a substantial mass renormalization factor $m^{*}/m_{\rm b}$
for all frequencies, ranging from 1.4 to 4.5, signifying strong electron
correlations. We note that $m^{*}/m_{\rm b}$ is larger for the hole pocket
($F_{3}$) compared to the electron pockets with linear band-crossings,
possibly reflecting the different impact of correlation effects on the
conventional and Dirac-like bands.

The inclusion of $U$ up to 3 eV in our calculations leads to a considerable
modification in the Fermi surface topology, as shown in Supplementary Fig. 4.
The best agreement with experiment, however, is found when $U$ = 0, indicating
that the influence of electron correlations cannot be effectively captured for
complex iridates using the GGA+SOC+$U$ approach [23]. A more accurate
accounting of correlation effects will likely require the implementation of
computationally demanding dynamical mean-field theory calculations, which will
be the subject of a future investigation. Nevertheless, the remarkable
agreement between the LT model and experiment demonstrates that correlation
effects play only a minor role in determining the Fermi surface topology, as
found for related systems such as $o$-SIO3 [23] and the topological semimetal
ZrSiS [29].

The anomalous nature of the electronic state in $m$-SIO3 is revealed through
its transport and thermodynamic properties at low temperatures. $\rho_{xx}(T)$
below 15 K follows an unusual $T+T^{2}$-dependence, as revealed by the finite
intercept in d$\rho_{xx}$/d$T$ at $T=0$ (Fig. 4a, b). The $T$-linear term
$A_{1}$ = 0.1 $\mu\Omega$ cm K-1, though small, persists down to the lowest
temperatures, possibly indicating the presence of low-energy critical
fluctuations. By contrast, the specific heat and magnetic susceptibility of
$m$-SIO3, shown in Fig. 4c, d, are more characteristic of a paramagnetic
Fermi-liquid state. Using the coefficient of the $T^{2}$ term $A_{2}$ = 0.048
$\mu\Omega$ cm K-2 and the electronic specific heat coefficient $\gamma_{0}$ =
4.47 (7) mJ mol-1 K-2, the corresponding Kadowaki-Woods ratio ($R_{\rm KW}$)
is found to be $A_{2}/\gamma_{0}^{2}\approx$ 2400 $\mu\Omega$ cm K2 mol2 J-2
or
$A_{2}/\gamma_{V}^{2}\approx$ 3.5 $\mu\Omega$ cm K2 cm6 J-2, where
$\gamma_{V}$ is the volume form of $\gamma_{0}$ [30]. These values, though
high relative to other correlated metals [30, 31, 32], are nonetheless
consistent with theoretical estimates derived for a 3D Fermi surface [30] (see
Supplementary Discussion for details), implying that electron-electron
scattering dominates the transport behaviour at low $T$. Recently, it has been
shown that the strength of electron interactions correlates with $R_{\rm KW}$
in the most general case [32], and overall consistency in our analysis
indicates that electron-electron scattering is responsible for the large
$R_{\rm KW}$ and the highly renormalized $m^{*}$ in $m$-SIO3, possibly due to
the SOC-renormalized bandwidth and/or proximity to a quantum critical point
[33, 34].

The topological character of the electronic bands in $m$-SIO3 is examined by
extracting the phase factor of the SdH oscillations (see Supplementary Fig.
5). While we find a non-trivial phase factor $\approx\rm{\pi}$, as expected
for a topologically non-trivial band, recent theoretical findings have
suggested that the phase factor in quantum oscillations can comprise of
multiple contributions and therefore cannot be considered as a conclusive
evidence for a non-trivial band topology [35]. Nonetheless, we emphasise that
the topological non-trivial character for the electron band in $m$-SIO3 is
further supported by the good agreement between the DFT calculations – with
their topologically-protected Dirac points at M and A – and the observed
quantum oscillation angular dependence.

Recently, transport signatures of Dirac quasiparticles have been reported in a
closely-related iridate CaIrO3 [36, 37]. The resultant Fermi surface areas and
effective masses, however, were found to be small ($F_{\rm SdH}$ = 3.2 $-$
11.2 T and $m^{*}$ = 0.12 $-$ 0.31 $m_{\rm e}$), suggesting that the Fermi
level in CaIrO3 is much closer to the Dirac point and that electron
correlations are considerably weaker than in $m$-SIO3. A recent photoemission
experiment on $o$-SIO3 thin films, on the other hand, found an electron-like
band with linear dispersion and parabolic hole-like bands with heavy
quasiparticle masses $m^{*}=2.4-6.0$ $m_{\rm e}$ [23], largely in agreement
with our findings on $m$-SIO3. Collectively, these findings identify SIO3 as a
rare example of a topological semimetal with enhanced electron correlations,
possibly induced through proximity to a Mott instability [9]. The extreme
sensitivity of its low-energy electronic structure to atomic displacement,
highlighted by the DFT calculations, further identifies $m$-SIO3 as a viable
platform on which to explore the Mott transition in a topological semimetal,
by tuning the relative strength of the electronic bandwidth and Coulomb
repulsion via doping or via structural tuning parameters such as hydrostatic
pressure or uniaxial strain.

## Methods

### Crystal growth and characterization

Single crystals of $m$-SIO3 were grown using the flux method described in Ref.
[38]. Low-temperature structural refinement using X-ray diffraction (XRD) was
performed at 13 K (see Supplementary Methods for details).

### Measurements

Electrical resistivity measurements were performed on Hall bars of ($600\times
100\times 60$) $\mu$m3 in size, cut from pristine crystals, with an ac current
of 3 mA applied along the bar identified as the crystallographic axis [100] by
Laue XRD. Magnetic fields up to 35 T were generated using a Bitter magnet at
the High Field Magnet Laboratory (HFML) in Nijmegen, the Netherlands. Zero-
field resistivity measurements were performed using a Cryogenic Free
Measurement System by Cryogenic Limited; specific heat and magnetic
susceptibility measurements were performed using a Physical Properties
Measurement System and a Magnetic Properties Measurement System XL,
respectively, by Quantum Design Inc.

### Band structure calculations

DFT band structure calculations were carried out using ab-initio VASP code
[39] with energy cut-off set to 400 eV and Brillouin zone sampled by the
$9\times 9\times 5$ Monkhorst-Pack $k$-point mesh. We used WANNIER90 package
[40] to construct the corresponding tight-binding (TB) Hamiltonian and ensured
that the resulting VASP and TB electronic band structures are identical near
the Fermi level. The visualization of Fermi surfaces was done using XCrysDen
code [41] and the oscillation frequency expected from the DFT Fermi surface
model was calculated using the SKEAF code [42].

## Acknowledgements

We gratefully acknowledge useful discussions with A. Rost and D. F. McMorrow.
We would also like to thank G. Stenning and D. Nye for help with the
instruments in the Materials Characterisation Laboratory at the ISIS Neutron
and Muon Source, Kuang-Yu Samuel Chang and Roos Leenen for technical
assistance with the DFT calculations, and Sebastian Bette for XRD
characterizations. We acknowledge the support of the HFML-Radboud University
(RU)/Netherlands Organisation for Scientific Research (NWO), a member of the
European Magnetic Field Laboratory. This work is part of the research program
Strange Metals (Grant 16METL01) of the former Foundation for Fundamental
Research on Matter, which is financially supported by the NWO and the European
Research Council (ERC) under the European Union s Horizon 2020 research and
innovation programme (Grant Agreement No. 835279-Catch-22). We gratefully
acknowledge support from the UK Engineering and Physical Sciences research
council, grant EP/N034694/1. We acknowledge collaborative support from A.S.
Gibbs, D. Fortes and the ISIS Crystallography Group for making available the
193Ir for the isotope work. Experiments at the ISIS Neutron and Muon Source
were supported by a beamtime allocation RB1990395,
DOI:10.5286/ISIS.E.RB1990395,? from the Science and Technology Facilities
Council. The work of D. P. and V. M. was supported by Act 211 Government of
the Russian Federation, contract 02.A03.21.0006.

## References

## References

  * [1] Keimer, B. & Moore, J. E.  The physics of quantum materials.  _Nat. Phys._ 13, 1045–1055 (2017).
  * [2] Dzero, M., Sun, K., Galitski, V. & Coleman, P.  Topological Kondo insulators.  _Phys. Rev. Lett._ 104, 106408 (2010).
  * [3] Pezzini, S. _et al._ Unconventional mass enhancement around the Dirac nodal loop in ZrSiS.  _Nat. Phys._ 14, 178–183 (2018).
  * [4] Kang, M. _et al._ Dirac fermions and flat bands in the ideal kagome metal FeSn.  _Nat. Mater._ 19, 163–169 (2020).
  * [5] Shi, W. _et al._ A charge-density-wave topological semimetal.  _Nat. Phys._ 17, 381–387 (2021).
  * [6] Tian, Z. _et al._ Field-induced quantum metal-insulator transition in the pyrochlore iridate Nd2Ir2O7.  _Nat. Phys._ 12, 134–138 (2016).
  * [7] Rau, J. G., Lee, E. K.-H. & Kee, H.-Y.  Spin-orbit physics giving rise to novel phases in correlated systems: iridates and related materials.  _Annu. Rev. Condens. Matter Phys._ 7, 195–221 (2016).
  * [8] Cao, G. & Schlottmann, P.  The challenge of spin-orbit-tuned ground states in iridates: a key issues review.  _Rep. Prog. Phys._ 81, 042502 (2018).
  * [9] Moon, S. J. _et al._ Dimensionality-controlled insulator-metal transition and correlated metallic state in $5d$ transition metal oxides Srn+1IrnO3n+1($n$ = 1, 2 and $\infty)$.  _Phys. Rev. Lett._ 101, 226402 (2008).
  * [10] Kim, B. J. _et al._ Novel J${}_{\rm eff}=1/2$ Mott state induced by relativistic spin-orbit coupling in Sr2IrO4.  _Phys. Rev. Lett._ 101, 076402 (2008).
  * [11] Kim, Y. K. _et al._ Fermi arcs in a doped pseudospin-1/2 Heisenberg antiferromagnet.  _Science_ 345, 187–190 (2014).
  * [12] Kim, Y. K., Sung, N. H., Denlinger, J. D. & Kim, B. J.  Observation of a $d$-wave gap in electron-doepd Sr2IrO4.  _Nat. Phys._ 12, 37–41 (2015).
  * [13] de la Torre, A. _et al._ Collapse of the Mott gap and emergence of a nodal liquid in lightly doped Sr2IrO4.  _Phys. Rev. Lett._ 115, 176402 (2015).
  * [14] Battisti, I. _et al._ Universality of pseudogap and emergent order in lightly doped Mott insulators.  _Nat. Phys._ 13, 21–25 (2017).
  * [15] Zhao, L. _et al._ Evidence of an odd-parity hidden order in a spin-orbit coupled correlated iridate.  _Nat. Phys._ 12, 32–36 (2015).
  * [16] Ding, Y. _et al._ Pressure-induced confined metal from the Mott insulator Sr3Ir2O7.  _Phys. Rev. Lett._ 116, 216402 (2016).
  * [17] Wu, F.-X. _et al._ Metal-insulator transition in SrIrO3 with strong spin-orbit interaction.  _J. Phys.: Condens. Matter_ 25, 125604 (2013).
  * [18] Biswas, A., Kim, K.-S. & Jeong, Y. H.  Metal insulator transition in perovskite SrIrO3 thin films.  _J. Appl. Phys._ 116, 213704 (2014).
  * [19] Groenendijk, D. J. _et al._ Spin-orbit semimetal SrIrO3 in the two-dimensional limit.  _Phys. Rev. Lett._ 119, 256403 (2017).
  * [20] Carter, J.-M., Shankar, V. V., Zeb, M. A. & Kee, H.-Y.  Semimetal and topological insulator in perovskite iridates.  _Phys. Rev. B_ 85, 115105 (2012).
  * [21] Zeb, M. A. & Kee, H.-Y.  Interplay between spin-orbit coupling and Hubbard interaction in SrIrO3 and related $Pbnm$ perovskite oxides.  _Phys. Rev. B_ 86, 085149 (2012).
  * [22] Takayama, T., Yaresko, A. N. & Takagi, H.  Monoclinic SrIrO3 – a Dirac semimetal produced by non-symmorphic symmetry and spin-orbit coupling.  _J. Phys.: Condens. Matter_ 31, 074001 (2019).
  * [23] Nie, Y. F. _et al._ Interplay of spin-orbit interactions, dimensionality, and octahedral rotations in semimetallic SrIrO3.  _Phys. Rev. Lett._ 114, 016401 (2015).
  * [24] Manca, N. _et al._ Balanced electron-hole transport in spin-orbit semimetal $sriro_{3}$ heterostructures.  _Phys. Rev. B_ 97, 081105(R) (2018).
  * [25] Sen, K. _et al._ Strange semimetal dynamics in SrIrO3.  _Nat. Commun._ 11, 4270 (2020).
  * [26] Shoenberg, D.  _Magnetic oscillations in metals_ (Cambridge University Press, 1984).
  * [27] Qasim, I., Kennedy, B. J. & Avdeev, M.  Synthesis, structures and properties of transition metal doped SrIrO3.  _J. Mater. Chem. A_ 1, 3127–3132 (2013).
  * [28] van Thiel, T. C. _et al._ Coupling lattice instabilities across the interface in ultrathin oxide heterostructures.  _ACS Mater. Lett._ 2, 389–394 (2020).
  * [29] Mueller, C. S. A. _et al._ Determination of the Fermi surface and field-induced quasiparticle tunneling around the Dirac nodal loop in ZrSiS.  _Phys. Rev. Res._ 2, 023217 (2020).
  * [30] Hussey, N. E.  Non-generality of the Kadowaki Woods ratio in correlated oxides.  _J. Phys. Soc. Jpn._ 74, 1107–1110 (2005).
  * [31] Jacko, A. C., Fjærestad, J. O. & Powell, B. J.  A unified explanation of the Kadowaki Woods ratio in strongly correlated metals.  _Nat. Phys._ 5, 422–425 (2009).
  * [32] Cavanagh, D. C., Jacko, A. C. & Powell, B. J.  Breakdown of the universality of the Kadowaki-Woods Ratio in multi-band metals.  _Phys. Rev. B_ 92, 195138 (2015).
  * [33] Cao, G. _et al._ Non-Fermi-liquid behavior in nearly ferromagnetic SrIrO3 single crystals.  _Phys. Rev. B_ 76, 100402(R) (2007).
  * [34] Groenendijk, D. J. _et al._ Anisotropic magnetoresistance in spin-orbit semimetal SrIrO3.  _Eur. Phys. J. Plus_ 135, 627 (2020).
  * [35] Alexandradinata, A., Wang, C., Duan, W. & Glazman, L.  Revealing the topology of fermi-surface wave functions from magnetic quantum oscillations.  _Phys. Rev. X_ 8, 011027 (2018).
  * [36] Fujioka, J. _et al._ Strong-correlation induced high-mobility electrons in Dirac semimetal of perovskite oxide.  _Nat. Commun._ 10, 362 (2019).
  * [37] Yamada, R. _et al._ Large variation of Dirac semimetal state in perovskite CaIrO3 with pressure-tuning of electron correlation.  _Phys. Rev. Lett._ 123, 216601 (2019).
  * [38] Hunter, E. C.  _An exploration of novel correlated electronic states in 5d transition metal oxides_.  Ph.D. thesis, The University of Edinburgh (2016).
  * [39] Kresse, G. & Furthmuller, J.  Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.  _Phys. Rev. B_ 54, 11169 (1996).
  * [40] Mostofi, A. A. _et al._ An updated version of wannier90: A tool for obtaining maximally-localised Wannier functions.  _Comput. Phys. Commun._ 185, 2309–2310 (2014).
  * [41] Kokalj, A.  XCrySDen – a new program for displaying crystalline structures and electron densities.  _J. Mol. Graph. Model._ 17, 176–179 (1999).
  * [42] Rourke, P. M. C. & Julian, S. R.  Numerical extraction of de Haas-van Alphen frequencies from calculated band energies.  _Comput. Phys. Commun._ 183, 324 (2012).