# Google COVID-19 Community Mobility Reports: Anonymization Process
Description (version 1.1)

Ahmet Aktay<EMAIL_ADDRESS>Shailesh Bavadekar  community-
<EMAIL_ADDRESS>Gwen Cossoul<EMAIL_ADDRESS>John Davis
<EMAIL_ADDRESS>Damien Desfontaines  community-
<EMAIL_ADDRESS>Alex Fabrikant<EMAIL_ADDRESS>Evgeniy
Gabrilovich<EMAIL_ADDRESS>Krishna Gadepalli  community-
<EMAIL_ADDRESS>Bryant Gipson<EMAIL_ADDRESS>Miguel
Guevara<EMAIL_ADDRESS>Chaitanya Kamath  community-
<EMAIL_ADDRESS>Mansi Kansal<EMAIL_ADDRESS>Ali Lange
<EMAIL_ADDRESS>Chinmoy Mandayam<EMAIL_ADDRESS>Andrew Oplinger<EMAIL_ADDRESS>Christopher Pluntke
<EMAIL_ADDRESS>Thomas Roessler<EMAIL_ADDRESS>Arran Schlosberg<EMAIL_ADDRESS>Tomer Shekel  community-
<EMAIL_ADDRESS>Swapnil Vispute<EMAIL_ADDRESS>Mia Vu
<EMAIL_ADDRESS>Gregory Wellenius  community-
<EMAIL_ADDRESS>Brian Williams<EMAIL_ADDRESS>Royce J
Wilson<EMAIL_ADDRESS>
###### Abstract

This document describes the aggregation and anonymization process applied to
the initial version of Google COVID-19 Community Mobility Reports (published
at http://google.com/covid19/mobility on April 2, 2020), a publicly available
resource intended to help public health authorities understand what has
changed in response to work-from-home, shelter-in-place, and other recommended
policies aimed at flattening the curve of the COVID-19 pandemic. Our
anonymization process is designed to ensure that no personal data, including
an individual’s location, movement, or contacts, can be derived from the
resulting metrics.

The high-level description of the procedure is as follows: we first generate a
set of anonymized metrics from the data of Google users who opted in to
Location History. Then, we compute percentage changes of these metrics from a
baseline based on the historical part of the anonymized metrics. We then
discard a subset which does not meet our bar for statistical reliability, and
release the rest publicly in a format that compares the result to the private
baseline.

COVID-19 Community Mobility Reports provide insights into changes in mobility
patterns. These reports use anonymized, aggregated data to chart movement
trends over time by geography, as well as by place categories, showing trends
over several weeks. This works in a similar way to existing Google products
and features. For example, Google Maps uses aggregated, anonymized data to
show how busy certain types of places are, including when a local business
tends to be the most crowded. Public health officials have suggested this same
type of aggregated, anonymized data could also be helpful as they make
critical decisions to combat COVID-19.

The COVID-19 Community Mobility Reports provide insights into what has changed
in response to work-from-home, stay-at-home, and other recommended policies
aimed at flattening the curve of the COVID-19 pandemic. They analyze trends in
visits made to high-level categories of places, including workplaces, retail
and recreational venues, groceries and pharmacies, parks, transit centers, and
places of residence. Each version of the report will show trends over several
weeks, with the most recent data representing 48 hours prior.

As explained in greater technical detail below, the anonymization process for
these reports includes differential privacy [1], which is well-suited to
produce analytics in contexts where the categories of data are known in
advance. Our rigorous approach intentionally adds random noise to metrics in a
way that maintains both users’ privacy and the overall accuracy of the
aggregated data.

This paper is structured as follows: we introduce our method to produce
anonymized metrics with differential privacy. We then explain how we post-
process the anonymized metrics to generate the reports. Figure 1 summarizes
the anonymization process.

OriginaldataPer-user contributionfor non-baseline daysPer-user contributionfor
baseline daysDaily metricsfor non-baseline daysDaily metricsfor baseline
daysNoisy metricsfor non-baseline daysNoisy metricsfor baseline daysPercentage
changesbetween daily metricsand baseline metricsPer-weekday medianof baseline
metricsPublished data extract and bound extract and bound aggregate across
users aggregate across users add noise add noise compute median per weekday
privacy boundary calculate ratio filter unreliable metrics Figure 1: System
diagram of the metrics computation and anonymization process

##  1 Definitions

#### Location History users

The metrics in these reports are based on the data of Google users who have
opted in to Location History [2], (“LH users”), a feature which is off by
default.

#### Differential Privacy [3]

Let $\varepsilon$ be a positive real number and A be a randomized algorithm
that computes a metric. In the context of this report, A is considered
$\varepsilon$-differentially private if for all input datasets $D_{1}$ and
$D_{2}$ such that $D_{2}$ can be obtained from $D_{1}$ by adding or removing a
single user’s data in a single day, and for all subsets of $S\in\text{im}{A}$:

$Pr[A(D_{1})\in S]\leq e^{\varepsilon}\cdot Pr[A(D_{2})\in S].$

#### Granularity levels

The metrics are aggregated per day and per geographic area. There are three
levels of geographic areas; in this paper, we call these granularity levels.

  * •

Granularity level 0 corresponds to metrics aggregated by country / region.

  * •

Granularity level 1 corresponds to metrics aggregated by top-level
geopolitical subdivisions (e.g. US states).

  * •

Granularity level 2 corresponds to metrics aggregated by higher-resolution
granularity (e.g. U.S. Counties).

Granularity levels 1 and 2 are defined differently in different countries, to
account for knowledge of local public-health needs. Note that in general, the
geographic area represented gets smaller as the granularity level increases.
No metrics are published for geographic regions smaller than 3km2.

##  2 Generating anonymized metrics

We are releasing aggregated, anonymized data that is designed to ensure that
no personal data, including an individual’s location, movement, or contacts,
can be derived from the resulting metrics. To that end, we anonymize the
statistics with differential privacy. We query the underlying data using our
open-source differential privacy library [4], which adds Laplace noise [5] to
protect each metric with differential privacy.

###  2.1 Daily visits in public places

We count the number of unique LH users who visited a public place of a given
category in a given day at each granularity level. There are seven different
categories derived from the data: retail, recreation, eateries (reported as
part of “Retail & recreation”); groceries, pharmacies; transit; and parks. We
add Laplace noise to each count according to the following table.

Granularity level | Scale of Laplace noise |  Corresponding $\varepsilon$ parameter
---|---|---
0 | $1/0.11\approx 9.09\quad(\sigma\approx 12.86)$ | 0.11
1 | $1/0.11\approx 9.09\quad(\sigma\approx 12.86)$ | 0.11
2 | $1/0.22\approx 4.55\quad(\sigma\approx 6.43)$ | 0.22
Table 1:  Noise parameters used for the daily visits in public places metrics

For each location (at all geographic levels), each LH user can contribute at
most once to each category. We also bound the contribution of each LH user to
4 ⟨category,location⟩ pairs per day and per geographic level, using a process
similar to the one described in this paper [6]: if an LH user contributes to
more than 4 pairs in a given day and given geographic level, we randomly
select 4 of them, and discard the others.

For example, suppose that on the same day, an LH user goes to public places in
all 7 categories in two distinct neighboring countries. This makes a total of
14 ⟨category,location⟩ pairs at country level. We would randomly discard 10 of
these pairs when computing country-level statistics.

This process does not significantly affect data accuracy: in the US, at county
level, $99\%$ of LH users contribute 3 or fewer ⟨category,place⟩ pairs per day
on average. Thus, each daily place visit is protected by differential privacy
with $\varepsilon=0.44$. These multiple metrics apply to the same dataset, so
standard composition results apply, and the total daily contribution of each
user is protected by differential privacy with a maximum of
$\varepsilon=1.76$.

###  2.2 Residential

For the purposes of this analysis, we use signals like relative frequency,
time and duration of visits to calculate metrics related to places of
residence. We calculate an average amount of time spent at places of residence
for LH users in hours. This computation is performed for each day and
geographic area, using the same algorithm as the differentially private mean
mechanism from our open-source library [7]. This mechanism works as follows:

  * •

We compute the amount of time spent at place of residence in a given day and
geographic area in hours by summing up the individual values per user offset
by 12, so all individual values fall into the range $[-12;12]$. We then add
Laplace noise to this sum; the scale of the noise is indicated in the table
below. We denote the real sum $s$, and noisy sum $s_{n}$.

  * •

We compute the count of unique users who spent any time at residences in a
given day and geographic area. We then add Laplace noise to this count; the
scale of the noise is indicated in the table below. We refer to the real count
$c$, and the noisy count $c_{n}$.

  * •

Finally, we compute the ratio $s_{n}/c_{n}$ for each day and each geographic
area, add 12 as offset, and clamp it to the range $[0,24]$ hours/day.

For example, at county-level, $s_{n}$ is obtained by first sampling a random
number from a Laplace distribution of scale 109.1, and then adding that number
to $s$. In the table below, we also indicate the standard deviation $\sigma$
of the noise added to each value.

Granularity level | Scale of Laplace noise: | Scale of Laplace noise: | Corresponding
---|---|---|---
| sum (total hours/day) | count (number of users) |  $\varepsilon$ parameter
0 | $12/0.055\approx 218.2\quad(\sigma\approx 308.6)$ | $1/0.055\approx 18.2\quad(\sigma\approx 25.71)$ | 0.11
1 | $12/0.055\approx 218.2\quad(\sigma\approx 308.6)$ | $1/0.055\approx 18.2\quad(\sigma\approx 25.71)$ | 0.11
2 | $12/0.110\approx 109.1\quad(\sigma\approx 154.3)$ | $1/0.110\approx 9.09\quad(\sigma\approx 12.86)$ | 0.22
Table 2:  Noise parameters used for the residential metrics

Each user can contribute to at most one region per granularity level, which
protects these metrics by differential privacy with $\varepsilon=0.44$ total
budget across all granularities. A description of the differentially private
mean mechanism implemented and a proof of its privacy guarantees is described
in [8] (Algorithm 2.4).

###  2.3 Workplaces

For the purposes of this analysis, we use signals like relative frequency,
time and duration of visits to calculate metrics related to places of
residence and places of work of LH users. We calculate how many LH users spent
more than 1 hour at their places of work. This computation is performed for
each day and geographic area. Then, we add Laplace noise to each count
according to the following table.

Geographic level | Scale of Laplace noise |  Corresponding $\varepsilon$ parameter
---|---|---
0 | $1/0.11\approx 9.09\quad(\sigma\approx 12.86)$ | 0.11
1 | $1/0.11\approx 9.09\quad(\sigma\approx 12.86)$ | 0.11
2 | $1/0.22\approx 4.55\quad(\sigma\approx 6.43)$ | 0.22
Table 3:  Noise parameters used for the work places metrics

The count is aggregated by places of residence of LH users. Since each user
can contribute to at most one geographic area per granularity level, these
metrics are protected by differential privacy with $\varepsilon=0.44$.

##  3 Generating the report from the anonymized metrics

The metrics described above are generated for each day, starting on
2020-01-01. They are then used to generate the percentage changes relative to
day of the week published in the reports. All operations described below use
only the output of the differentially private mechanisms described in the
previous section; so they do not consume any privacy budget.

#### Additional privacy protections

We discard all metrics for which the geographic region is smaller than 3km2,
or for which the differentially private count of contributing users (after
noise addition) is smaller than 100. Geographic regions smaller than 3km2 may
be merged such that the union of their area is above the 3km2 threshold. This
merging does not occur across country boundaries, except for the Vatican City
and Italy.

###  3.1 Computing percentage changes from a baseline

For each individual metric generated using the mechanisms described above, we
compute the ratio between the metric for a given day D and the same metric
computed for the baseline period. The reference baseline is defined in the
following way.

  * •

We consider the 5-week range from 2020-01-03 through 2020-02-06. This ranged
is fixed.

  * •

Within this 5-week range, we consider the 5 days with the same day of week as
$d$. For example, if $d$ is 2020-03-20, $d$ is a Friday, so we consider the 5
Fridays in this 5-week range (Jan 3 to Jan 31, inclusive).

  * •

We compute the median of the differentially private metrics for these 5
baseline days.

  * •

This median metric is the baseline metric for $d$.

We then compute and publish the ratio between the metric for $d$ and the
baseline metric, as a percentage.

###  3.2 Removing unreliable metrics

In some regions, the noise added to obtain differential privacy can reduce the
confidence that we are capturing a meaningful change, typically when there is
not a lot of data for the metric. When, because of this uncertainty, the
percentage change for one of these metrics has a $5\%$ chance (or higher) of
being wrong by more than $\pm 10$ absolute percentage points, we do not
publish it and instead include an asterisk denoting that there is not enough
data available to present privacy-safe information. More precisely:

  * •

Before releasing a ratio metric/baseline, we compute $97.5\%$ confidence
intervals for the metric and its baseline. Let us denote [$m_{\min}$,
$m_{\max}$] and [$b_{\min}$, $b_{\max}$] these respective confidence
intervals.

  * •

We compute the ratios $m_{\min}/b_{\max}$ and $m_{\max}/b_{\min}$.

  * •

If one of these ratios differs from the differentially private ratio by more
than 10 absolute percentage points, we do not publish the corresponding
percentage changes.

If the last condition is not satisfied, then the probability of being wrong by
more than 10 absolute percentage points in each direction is lower than 2.5%.
By union bound, this means that there is at most a 5% risk of being wrong by
more than 10 absolute percentage points. Note that the confidence intervals
are based on an already differentially private value and on public data (the
scale and shape of the noise), so no privacy budget is consumed by this
operation.

##  4 Note on $\delta$

We are generating a fixed set of metrics (all possible combinations of
geographic regions, days within the periods, and public place categories), and
we are also adding noise to zero-valued metrics. As such, the process outlined
above is $\varepsilon$-differentially private with $\delta=0$.

##  5 Improving the accuracy of metrics over time

We are continuously making improvements to the underlying computation of the
metrics to improve their accuracy over time. These updates can introduce a
shift in their values, which can skew comparisons over time compared to the
baseline period if they are not accounted for, since those improvements are
not applied to the baseline period (to avoid republishing the data). To roll
out these changes with minimal impact on the overall privacy budget, we use
_scaling factors_. Rather than recomputing a metric directly, which would
require a larger value of $\varepsilon$, we use the following process.

  1. 1.

We define groups of metrics, on which the effect of the update is uniform. For
example, for a specific metric and location, we could consider all days in a
given period, or group days by weekday within this period (e.g. Parks metric
across all Tuesdays in June in the US, or Workplaces metric across week-ends
in August in each region of granularity level 1).

  2. 2.

For each group, we sum the noisy metrics already generated with the previous
computation logic. Call this sum $s_{g}$.

  3. 3.

For each group, we sum the metrics recomputed for the same given period with
the new computation logic, then we add noise to it with a smaller privacy
budget (typically, 10%), proportional to the budget we use for the
corresponding region granularity. Call this noisy sum $s_{n}$.

  4. 4.

For dates following the given period, we multiply the metrics generated with
the new logic by the scaling factor $s_{n}/s_{g}$ (if we are scaling the
baseline) or by its inverse $s_{g}/s_{n}$ (if we are scaling the daily
counts).

Grouping multiple metrics together to compute this scaling factor is the key
insight that allows us to use a much smaller privacy budget for step 3, and
reuse already generated metrics in step 2.

## References

  * [1] Enabling developers and organizations to use differential privacy.  https://developers.googleblog.com/2019/09/enabling-developers-and-organizations.html.
  * [2] Location History.  https://support.google.com/accounts/answer/3118687?hl=en.
  * [3] Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith.  Calibrating Noise to Sensitivity in Private Data Analysis.  In TCC, volume 3876 of Lecture Notes in Computer Science, pages 265–284. Springer, 2006.
  * [4] Google’s C++ Differential Privacy Library.  https://github.com/google/differential-privacy.
  * [5] Laplace Distribution.  https://en.wikipedia.org/wiki/Laplace_distribution.
  * [6] Royce Wilson, Celia Yuxin Zhang, William Lam, Damien Desfontaines, Daniel Simmons-Marengo, and Bryant Gipson.  Differentially Private SQL with Bounded User Contribution.  In Proceedings on Privacy Enhancing Technologies Symposium (to appear), 2020.
  * [7] DP Library Bounded Mean.  https://github.com/google/differential-privacy/blob/master/differential_privacy/algorithms/bounded-mean.h.
  * [8] N. Li, M. Lyu, D. Su, and W. Yang.  Differential Privacy: From Theory to Practice.  Morgan & Claypool Publishers, 2016.


# What is the Standard Cosmological Model?

Eric V. Linder1,2 1Berkeley Center for Cosmological Physics & Berkeley Lab,
University of California, Berkeley, CA 94720, USA
2Energetic Cosmos Laboratory, Nazarbayev University, Nur-Sultan 010000,
Kazakhstan

###### Abstract

Reports of “cosmology in crisis” are in vogue, but as Mark Twain said, “the
report of my death was an exaggeration”. We explore what we might actually
mean by the standard cosmological model, how tensions – or their apparent
resolutions – might arise from too narrow a view, and why looking at the big
picture is so essential. This is based on the seminar “All Cosmology, All the
Time”.

##  I Introduction

The origin of the word “crisis” comes from “decision”, implying a framework –
a standard model – and paths forward that are being decided between. So to
assess whether some crisis exists and how severe it might be, we must first
explore what we mean by the standard cosmological model.

We discuss various definitions in Section II, then investigate a couple of oft
claimed tensions and their (un)resolutions in Section III. The primacy of
using all robust observations, at all redshifts – All Cosmology, All The Time
– is motivated in Section IV, and the concluding discussion is in Section V.

##  II What is Standard Cosmology?

What is the standard cosmological model depends very much on where one draws
the line on what is cosmology or the universe.

###  II.1 Level 1: Global Properties

One could define Level 1 as saying cosmology is the global properties of the
universe:

  * •

Connected – a signal can get from there to here and from then to now. There is
no discreteness (at least at this level).

  * •

Metric – we can figure out how far it is from there to here and from then to
now.

This gives the foundation and few would dispute that they are part of the
standard cosmological model; violation of either on the scale of the
observable universe would be revolutionary (and pretty exciting on smaller
scales as well!).

A third property comes from a plethora of observations, and is also generally
accepted as part of our standard cosmological model:

  * •

Homogeneity and Isotropy – on observable universe scales.

This then implies the Robertson-Walker metric and we are on familiar
territory! In particular, the Robertson-Walker metric has two characteristics
(independent of the theory of gravity):

  * •

Evolution – a scale factor $a(t)$.

  * •

Spatial curvature – a constant $k$ proportional to the Gaussian curvature of
space. It is worth emphasizing that a Robertson-Walker spacetime can have
spatial curvature but no spacetime curvature (though this does not describe
our universe), and spacetime curvature but no spatial curvature (which fits
observations pretty well).

Unlike, say, discreteness, we don’t have to go down to laboratory scales to
test for a breakdown in homogeneity and isotropy, and hence the Robertson-
Walker metric. However, one can show that these smaller scale deviations do
not significantly affect cosmological scales. This is worth emphasizing: for
the expansion and curvature, and for observations involving propagation of
light rays along a line of sight, and for growth within three dimensional
volumes, the effect from small scales generally is calculable and found to be
small [1, 2, 3, 4].

###  II.2 Level 2: History of the Universe

If one delves into the time structure, one could say that a Level 2 answer to
what is the standard cosmological model is to say cosmology is the history of
the universe. From observations we could list five key stages:

  * •

Early hot dense state – popularly, if nebulously, known as the “Big Bang”.
Here we consider it as a generic description and source of initial conditions,
and whether it occurs at the Planck energy, $10^{15}$ GeV, or $10^{3}$ GeV is
a (fascinating) detail, not a source of crisis.

  * •

Matter/antimatter asymmetry – I have absolutely nothing to say! While we know
the basic ingredients to deliver this (if not its magnitude) [5] one could
well argue this is a greater crisis than any other raised, and yet it is
rarely mentioned.

  * •

Radiation dominated era – primordial nucleosynthesis, degrees of freedom
$g_{\star}$ (neutrino decoupling, electron/positron annihilation), CMB
thermalization. Lots of good stuff!

  * •

Matter dominated era – CMB scattering, growth of structure (us!).

  * •

Cosmic acceleration – “dark energy”, fate of the universe?

As an aside, one of the most fascinating aspects of this is that the argument
could be made that the study of dark energy actually grew out of investigation
of the radiation dominated era. The radiation era spans such a huge range of
e-folds of expansion (nominally some 55, compared to the matter era’s 7, or
dark energy’s 0.5), and yet we know little of the details – did it stay
radiation dominated the entire time? We tend to simply assume this but we have
accurate windows on only tiny slivers of this era, around primordial
nucleosynthesis [6, 7, 8] and toward the end, just before matter domination
and recombination [9, 10, 11]. In 1979, Robert Wagoner [12] highlighted this
by considering what freedom there was of changing the equation of state of the
dominant energy density. This program of exploring the equation of state at
various epochs was picked up by two of his students and in the 1980s was
applied to the late universe and became dark energy cosmology, developing both
the cosmological model and observational probes, first for particular values
of the equations of state [13, 14] and then for a general equation of state
[15, 16].

###  II.3 Level 3: Stuff In the Universe

A more specific view is that cosmology is the stuff in the universe. After
all, this is what is actually observed. This would include:

  * •

Cosmic microwave background radiation (CMB) – CMB structure (anisotropies,
polarization, spectral distortions) is a rich probe of both history (Level 2,
including initial conditions such as adiabatic perturbations) and the other
contents (Level 3, e.g. matter stuff through scattering, gravitational
potentials). And of course it provides strong evidence through isotropy (and
homogeneity through the CMB felt by distance objects) for Level 1 cosmology.

  * •

Large scale structure – The continuous fields of matter: the density field,
velocity field, acceleration (gravity) field – generally as probed by
individual sources. While these fields are related, the relations do test the
framework and each has particular incisive elements so they can be considered
distinct; plus, for cosmological purposes they are at very different stages of
observational development.

  * •

Other – As probes of the standard cosmological model, observations of other
stuff such as neutrinos, gravitational waves, exotica (e.g. topological
defects) have not yet reached the same stage of having a major impact, though
this would be an exciting development.

###  II.4 Turtles All the Way Down

Some researchers extend the standard model of cosmology further and further
down in detail, to “the stuff in the stuff in the universe”, e.g. aspects of
galaxies, galaxy clusters, generation of various particles and fields (e.g.
neutrinos, gravitational waves, etc.). Others will stretch to “the properties
of the stuff in the stuff in the universe”, e.g. cuspy cores of galaxies,
tidal streams, Cepheid pulsations, etc.

Where one draws the line between the standard cosmological model and all else
– call it astrophysics for simplicity – is a personal choice. But it is rare
that rain next Tuesday in some place one had not predicted it throws the
standard meteorological model into crisis. It is not impossible, but it is
useful to remember that some particular tension has to work its way from Level
$N$ to the fundamental foundations.

###  II.5 Cosmologing Is Hard

But…cosmologing is hard. The properties of the stuff in the stuff affect how
and what we learn about the more fundamental stuff. We are then faced with the
puzzle of whether we have fully understood the stuff in the stuff (or
mismeasured its properties) or whether indeed it propagates cleanly to
impacting the standard cosmological model. Let us take two examples.

Example 1. Suppose one measured the CMB temperature at scale factor $a$
(redshift $z=a^{-1}-1$) and found that

$T_{\rm CMB}(z)\neq T_{\rm CMB}(0)\times(1+z)\,?$ (1)

What should one conclude: that the universe is not adiabatically expanding, or
that there is some systematic error (e.g. molecular collisional excitations)?
How much effort, and what proportion of the literature, should be dedicated to
investigating systematics before declaring a crisis in the standard
cosmological model?

Example 2. Suppose one measured the distance to an object (or set of objects)
at redshift $z$ in terms of both its luminosity distance and angular diameter
distance and found that the reciprocity relation is broken,

$d_{L}(z)\neq d_{a}(z)\times(1+z)^{2}\,?$ (2)

If one wishes to conclude that this places the standard cosmological model in
crisis (rather than some systematic error in the data), one must give up some
foundational element, since this relation arises from a) metricity, b)
geodesic completeness, c) photons propagate on null geodesics, and d)
adiabatic expansion. (Note conservation of photon phase space density is a big
part of this, but not all.)

What level of systematics investigation should one (and the research
community) carry out before declaring an upending of the standard cosmological
model? To what extent should the proportion of systematics vs new physics
papers depend on the degree of new physics required? As the saying goes, if
you tell me you saw out the window a deer on the lawn, I might be willing to
consider how the deer got there, its effect on the shrubbery, etc., but if you
tell me you saw out the window a unicorn on the lawn, I might want more
evidence before devoting time to the puzzle.

We have in place solutions for how to handle conflicting expectations,
observations, and theories:

  * •

Rigorous data

  * •

Multiple, disparate probes

  * •

Crosschecks

  * •

Consistency at all cosmic times

  * •

Check the cosmic expansion history, cosmic growth history, and light
propagation (and soon gravitational wave propagation)

We explore how these might be applied to an example tension in the next
section.

##  III Past Tense, Present Tense?

A well known tension lies in current Hubble constant $H_{0}$ values deduced
from certain probes, taking all the data at face value. In addition to the
main puzzle, there are some beyond the surface:

  * •

Local measurements differ by some $\sim 2\sigma$ depending on method, i.e.
Cepheids vs tip of the red giant branch [17, 18].

  * •

The tension is emphatically not “early vs late” cosmology since baryon
acoustic oscillations (BAO) distances (together with primordial element
abundances [19, 20, 21] or marginalizing over or sidestepping the sound
horizon at the baryon drag epoch [22, 23]), i.e. without use of the primordial
CMB, gives the same answer as from the CMB.

  * •

Strong lensing time delays show a sharp transition between low and high
$H_{0}$ values around $z\sim 0.4$ [24, 25], albeit with a small sample.

While CMB data alone constrains $H_{0}$ tightly only within a $\Lambda$CDM
cosmological model, while allowing a considerable range of $H_{0}$ when the
dark energy equation of state $w$ differs from $-1$, it is extraordinarily
difficult from a combination of cosmic probes such as CMB+BAO or CMB+SN
(supernovae distances) to obtain $H_{0}>70$ [26] (we always write $H_{0}$ in
units of km/s/Mpc). Basically, $H_{0}>70$ requires a phantom dark energy
($w<-1$), which is disfavored by the above combinations of probes; see Figure
1.

Figure 1: 68.3% and 95.4% constraints on the $w_{0}$–$w_{a}$ plane in an $11$
parameter extended space, using Planck CMB data plus the R16 $H_{0}$ ($\sim
74$) prior. The top panel includes as well JLA supernova data, the middle
panel includes baryon acoustic oscillations data: both strongly prefer
$w_{0}\geq-1$. The bottom panel shows that $w_{0}\geq-1$ is most cleanly
achieved by shifting the $H_{0}$ prior to $H_{0}\leq 70$. (Adapted from [26].)

There are two basic loopholes one might try by changing the cosmic expansion
history – relax the tension in the present or the past.

  * •

Present tense (late time transition): arrange a very sharp phantom excursion
very close to the present so that higher redshift distances are not too
strongly affected.

  * •

Past tense (early time transition): arrange a lower sound horizon scale
$r_{\rm drag}$ with the Hubble parameter $H$ going up. Again, one must make it
a sharp transition – a spurt of extra early energy density to raise $H$, then
removing the early dark energy to preserve the agreement with CMB data.

We begin with the early time transition. The covariance between $r_{\rm drag}$
and $H_{0}$ has been known for a long time [27, 28, 29, 30, 11]. Using CMB
data, in 2013 Ref. [11] actually found evidence for an early time transition
and its effect on $H_{0}$! – see Figure 2. This has been resuscitated in many
many articles in the last couple years. However, early time transitions do not
really work in removing the Hubble constant tension (see, e.g., [31, 32, 33,
34]). On the one hand they cannot viably raise $H_{0}$ as far as the high
values favored by Cepheids, and on the other hand they generically violate
other aspects of the cosmological model as we discuss in Section IV.

Figure 2: Reconstruction of the expansion history deviations $\delta(a)=\delta
H^{2}/H^{2}_{\Lambda{\rm CDM}}$ from $\Lambda$CDM is shown, with the mean
value (solid line) and 68% uncertainty band (shaded area). Note the data
prefers early dark energy $\delta\neq 0$. The inset demonstrates that the
model independent reconstruction then prefers a higher $H_{0}$ than the
$\Lambda$CDM standard analysis. (Adapted from [11].)

The late time transition has the advantage that there is much less effect on
the CMB, and so more freedom for change. If one raises $H(z)$, distances will
change. To preserve distances (e.g. the distance to the CMB last scattering
surface, as well as BAO and SN distances), with a higher $H_{0}$ one needs a
smaller $H(z>0)$. This means less energy density. This could be from a smaller
matter density $\Omega_{m}$ (the matter density today as a fraction of the
critical density) but this is insufficient and one requires a smaller dark
energy density at some $z>0$. To give enough dark energy density today
(especially if the matter density is low and the total density is the critical
density) requires dark energy density to appear quite suddenly at low
redshifts. This is again the phantom regime $w<-1$. And again, such late time
phase transitions have been known for a long time (see, e.g., vacuum
metamorphosis [35, 36, 37]).

Since such late time transitions all have basically the same physical effect,
regardless of specific origin, let’s examine whether vacuum metamorphosis
removes the $H_{0}$ tension. Yes! but…. As the opening lines of the abstract
of [38] say, “We do obtain $H_{0}\approx 74$ km/s/Mpc from CMB+BAO+SN data in
our model, but that is not the point.” This – and essentially any late time
transition – model fails because there is more to the standard cosmological
model than $H_{0}$; it does not satisfy other probes. We discuss the details
in Section IV. For other possible issues with late time transitions see [31,
39, 40, 41].

##  IV All Cosmology, All The Time

So far we have considered only the expansion history. However, one must take
into account all cosmological probes, e.g. how deviations from a standard
cosmological model affect the growth of large scale structure.

In the vacuum metamorphosis case of the previous section, the combination of
probes CMB+BAO+SN produced $H_{0}\approx 74$. For a good fit to the CMB,
preserving $\Omega_{m}h^{2}$ means a low $\Omega_{m}\approx 0.27$. That can be
ok. However, it also gives a high amplitude for mass fluctuations,
$\sigma_{8}\approx 0.88$, which is quite high. This is due to the reduced dark
energy density needed to get the distances right, quite generally implying
greater matter domination and growth at higher redshifts. We can start to see
that we do indeed need “all cosmology, all the time” – use of all probes, over
the full cosmic history.

One might brush high $\sigma_{8}$ under the rug and say that with the lower
$\Omega_{m}$, one has $S_{8}\equiv\sigma_{8}(\Omega_{m}/0.3)^{0.5}\approx
0.83$, which might be workable for some probes, i.e. roughly as good as
$\Lambda$CDM. So we could say that vacuum metamorphosis gives $H_{0}\approx
74$ while not making any $S_{8}$ tension worse, as apparently seen in Figure
3.

However, $S_{8}$ and $H_{0}$ are focusing on a single time (the present) in
cosmic history. This is a Bad Idea. In the words of Lewis Carroll, “the rule
is, jam tomorrow and jam yesterday – but never jam today”. Figure 4 shows that
the apparent removal of the $H_{0}$ tension is moot, since Figure 3 is only a
tiny part of cosmic history, and when one does all cosmology, all the time –
taking into account both the cosmic expansion and cosmic growth over the span
of cosmic history as in Figure 4 – then generically late time transitions do
not work in giving a viable cosmological model.

Figure 3: Expansion and growth histories are here plotted simultaneously. At
$z=0$, hence these give $H_{0}$ and $f\sigma_{8}(0)$ ($\approx 0.52\,S_{8}$
for most viable cosmologies in general relativity; see [38]), shown by the
points. Excellent agreement with the Cepheid value of $H_{0}$ (green dashed
line), as opposed to the $\Lambda$CDM value (cyan dashed line) is obtained by
the flat vacuum metamorphosis (VM) late transition models with parameters fit
to CMB+BAO+SN, and the VM VEV model agrees well on $S_{8}$ as well. (Adapted
from [38].)  Figure 4: A very different situation occurs if one looks beyond
$H_{0}$ at the full conjoint evolutionary track of the various models, not
just the $z=0$ (lower) endpoints as shown in Figure 3 (which covers the small
range between the dashed lines). Over the histories the VM models diverge
considerably from $\Lambda$CDM. Curves extend from $z=0$ at the bottom to
$z=3$ at the top, and the points with error bars show the trajectory location
at $z=0.6$, 1, 2, where the error bars mimic 1% constraints on each axis
quantity to give a sense of separation between the curves. Solid curves are
for flat space, dotted curves include $\Omega_{k}$. (Adapted from [38].)

This is general because the late time transition, anchored by the present,
requires a lower dark energy density at higher redshifts and hence greater
growth. If you squeeze the model in one place, it will bulge out elsewhere and
fail to fit the array of probes. A similar general “no-go” reasoning holds for
early time transitions. The higher expansion rate damps the CMB perturbations,
so to preserve the fit to CMB data one requires a higher primordial curvature
perturbation amplitude, and this strengthens the later growth of the matter
perturbations, leading to high $\sigma_{8}$. (See also [42].)

##  V Discussion

Very generally, neither late time nor early time attempts to remove the
$H_{0}$ tension survive all cosmology, all the time. One needs to take into
account all the probes, at all redshifts. It is not just $H_{0}$, it is
$H(z)$. It is not just $\Omega_{m}$, it is $\Omega_{m}(z)$, and growth of
structure, light propagation, etc.

The standard cosmological model, whether one views it at Level 1, 2, or
whatever, is strong: I come to praise it, not to bury it.

But, arguendo, suppose we ignored the questions of Section II.5 and believed
all data blissfully free of systematics, so that the order (two orders?) of
magnitude difference in number of papers squeezing the standard cosmological
model vs investigating the data is right and proper. What could we do to
remove the tension? From the preceding section it is clear: we must break the
connection between the cosmic expansion history and growth history. We can do
this by a breakdown in general relativity, making gravity weaker so that
$\sigma_{8}$ becomes lower; we can do this by introducing new particle
interactions, again to suppress growth. Such changes will in turn alter other
cosmological probes, which must be considered, and we must be wary of a spiral
of epicycles. (See also [43].)

Basically we would need to break the standard expansion history to address the
$H_{0}$ tension, and need to break the standard growth history to keep the
consequences of the first break consistent with other probes. An interesting
possibility is to probe the connection of the growth of structure to the
cosmic expansion in a wholly new way. Gravitational wave standard sirens offer
this possibility, in part. A direct connection between siren distances and
matter growth was developed by [44].

The propagation of gravitational waves is an excellent probe of “spacetime
friction” – a combination of the Hubble friction from expansion and any time
variation of the gravitational strength; the growth of structure probes both
these quantities as well, in a different way. By comparing these two probes,
possibly plus light propagation which depends only on the expansion, we have
the potential to get a clear view of whether the connections between them in
the standard cosmological model hold – and if they do not, is there
consistency between the deviations in one probe and in another.

A bit more technically: while distances found through light propagation
$d_{EM}(z)=d_{L}(z)$ depend on the expansion $H(z)$, those determined through
gravitational wave propagation depend on both $H(z)$ and $\alpha_{M}(z)\equiv
d\ln M^{2}_{Pl}(a)/d\ln a$, the running of the Planck mass or inverse
gravitational strength. Thus, a measurement $d_{GW}(z)\neq d_{EM}(z)$ would
signal a deviation from general relativity (or systematics). Meanwhile, growth
of cosmic structure depends on $H(z)$ and the strength of gravity $G_{\rm
eff}(z)$. In theories where the gravitational strength is determined purely by
the Planck mass, $G_{\rm eff}(z)=1/M_{Pl}^{2}(z)$, then the circle is complete
and there is a tight relation between the three probes. (Some theories of
gravity have a further factor, the braiding that mixes the tensor and scalar
parts, loosening the relation; $G_{\rm eff}(z)$ can also affect light
propagation distances, and hence $H_{0}$ from strong lensing, giving a
redshift dependence to the derived $H_{0}$ [25].) Thus we can crosscheck
against systematics by seeing if the relation holds: a deviation in one probe
predicts a specific redshift dependence for a deviation in another probe.

The consistency check can be quantified with a new statistic combining the
probe measurements [45],

$D_{G}(a)\equiv\frac{\left[d_{L,GW}^{\rm MG}/d_{L}^{\rm
GR}\right](a)}{\left[f\sigma_{8}^{\rm MG}/f\sigma_{8}^{\rm GR}\right](a)}\,.$
(3)

For general relativity, this equals one for all redshifts. For a given
modified gravity (MG) model, it has a specific redshift dependence predicted.
Several examples are shown in Figure 5.

In summary, the standard cosmological model has extremely deep foundations and
multiple layers, and apparent surface blemishes may have very little to do
with the fundamental basis. Especially if tensions cannot be solved by one
“tooth fairy”, such as an early or late time transition, when confronted with
application of the principle of using all cosmology, all the time, but we are
instead led to a series of epicycles, we might recall the unicorn vs the deer
and pause.

Figure 5: The new $D_{G}$ statistic, using the complementarity of the
gravitational wave luminosity distance $d_{L,GW}$ and the cosmic matter growth
rate $f\sigma_{8}$, can clearly distinguish different classes of gravity. Each
class has a distinct shape in its redshift dependence $D_{G}(a)$, though the
curve amplitudes will scale with $G_{\rm eff}(z=0)$. General relativity has
constant $D_{G}=1$. (Adapted from [45].)

New probes and data covering more cosmic history will be essential in
exploring the standard cosmological model, at whatever level we define it, and
whether it will stand firm, need some patchwork, or be overturned. Each of the
levels has a wonderful array of areas for students to research, and make
significant and lasting contributions to what their students will learn as the
standard cosmological model.

###### Acknowledgements.

I thank the Asia-Pacific Center for Theoretical Physics for inviting me to
give this inaugural lecture in June 2020 for the series “Dark Energy in a Dark
Age”. This work is supported in part by the Energetic Cosmos Laboratory and by
the U.S. Department of Energy, Office of Science, Office of High Energy
Physics, under contract no. DE-AC02-05CH11231.

## References

  * [1] T. Futamase, An approximation scheme for constructing inhomogeneous universes in general relativity, MNRAS 237, 187 (1989)
  * [2] E.V. Linder, The linearly perturbed sky, MNRAS 243, 362 (1990)
  * [3] M.W. Jacobs, E.V. Linder, R.V. Wagoner, Obtaining the metric of our Universe, Phys. Rev. D 45, 3292 (1992)
  * [4] M.W. Jacobs, E.V. Linder, R.V. Wagoner, Green function for metric perturbations due to cosmological density fluctuations, Phys. Rev. D 48, 4623 (1993)
  * [5] A.D. Sakharov, Violation of CP Invariance, C asymmetry, and baryon asymmetry of the universe, Pisma Zh. Eksp. Teor. Fiz. 5, 32 (1967) [Sov. Phys. Usp. 34, 392 (1991)]
  * [6] S.M. Carroll, M. Kaplinghat, Testing the Friedmann equation: The expansion of the universe during big-bang nucleosynthesis, Phys. Rev. D 65, 063507 (2002) [arXiv:astro-ph/0108002]
  * [7] E. Masso, F. Rota, Primordial nucleosynthesis as a test of the Friedmann equation in the early Universe, Phys. Rev. D 68, 123504 (2003) [arXiv:astro-ph/0302554]
  * [8] S. Dutta, R.J. Scherrer, Big Bang nucleosynthesis with a stiff fluid, Phys. Rev. D 82, 083501 [arXiv:1006.4166]
  * [9] E.V. Linder, T.L. Smith, Dark Before Light: Testing the Cosmic Expansion History through the Cosmic Microwave Background, JCAP 1104, 001 (2011) [arXiv:1009.3500]
  * [10] J. Samsing, E.V. Linder, T.L. Smith, Model Independent Early Expansion History and Dark Energy, Phys. Rev. D 86, 123504 (2012) [arXiv:1208.4845]
  * [11] A. Hojjati, E.V. Linder, J. Samsing, New Constraints on the Early Expansion History, Phys. Rev. Lett. 111 041301 (2013) [arXiv:1304.3724]
  * [12] R.V. Wagoner, The Early Universe, in Proc. Les Houches Summer School 32, p. 395 (1980)
  * [13] M.S. Turner, Coherent scalar-field oscillations in an expanding universe, Phys. Rev. D 28, 1243 (1983)
  * [14] J.C. Charlton, M.S. Turner, Kinematic Tests of Exotic Flat Cosmological Models, ApJ 313, 495 (1987)
  * [15] E.V. Linder, Cosmological tests of generalized Friedmann models, Astron. Astroph. 206, 175 (1988)
  * [16] E.V. Linder, Light propagation in generalized Friedmann universes, Astron. Astroph. 206, 190 (1988)
  * [17] A.G. Riess, S. Casertano, W. Yuan, J.B. Bowers, L. Macri, J.C. Zinn, D. Scolnic, Cosmic Distances Calibrated to 1% Precision with Gaia EDR3 Parallaxes and Hubble Space Telescope Photometry of 75 Milky Way Cepheids Confirm Tension with $\Lambda$CDM, ApJL 908, L6 (2021) [arXiv:2012.08534]
  * [18] W.L. Freedman, B.F. Madore, T. Hoyt, I.S. Jang, R. Beaton, M.G. Lee, A. Monson, J. Neeley, J. Rich, Calibration of the Tip of the Red Giant Branch, ApJ 891, 57 (2020) [arXiv:2002.01550]
  * [19] G. Addison, D. Watts, C. Bennett, M. Halpern, G. Hinshaw, J. Weiland, Elucidating $\Lambda$CDM: Impact of Baryon Acoustic Oscillation Measurements on the Hubble Constant Discrepancy, ApJ 853, 119 (2018) [arXiv:1707.06547]
  * [20] A. Cuceu, J. Farr, P. Lemos, A. Font-Ribera, Baryon Acoustic Oscillations and the Hubble Constant: Past, Present and Future, JCAP 1910, 044 (2019) [arXiv:1906.11628]
  * [21] O.H.E. Philcox, M.M. Ivanov, M. Simonovíc, M. Zaldarriaga Combining Full-Shape and BAO Analyses of Galaxy Power Spectra: A 1.6% CMB-independent constraint on $H_{0}$, JCAP 2005, 032 (2020) [arXiv:2002.04035]
  * [22] E. Macaulay et al. (DES Collaboration), Mon. Not. Roy. Astron. Soc. 486, 2184 (2019) [arXiv:1811.02376]
  * [23] E.J. Baxter, B.D. Sherwin, Determining the Hubble Constant without the Sound Horizon Scale: Measurements from CMB Lensing, MNRAS 501, 1823 (2021) [arXiv:2007.04007]
  * [24] M. Millon et al., TDCOSMO. I. An exploration of systematic uncertainties in the inference of $H_{0}$ from time-delay cosmography, Astron. Astroph. 639, A101 (2020) [arXiv:1912.08027]
  * [25] K. Liao, A. Shafieloo, R.E. Keeley, E.V. Linder, Determining Model-independent $H_{0}$ and Consistency Tests, ApJL 895, L29 (2020) [arXiv:2002.10605]
  * [26] E. DiValentino, A. Melchiorri, E.V. Linder, J. Silk, Constraining Dark Energy Dynamics in Extended Parameter Space, Phys. Rev. D 96, 023523 (2017) [arXiv:1704.02153]
  * [27] G. Efstathiou, J. Bond, Cosmic Confusion: Degeneracies among Cosmological Parameters Derived from Measurements of Microwave Background Anisotropies, MNRAS 304, 75 (1999) [arXiv:astro-ph/9807103]
  * [28] D.J. Eisenstein, M.J. White, Theoretical uncertainty in baryon oscillations, Phys. Rev. D 70, 103523 (2004) [arXiv:astro-ph/0407539]
  * [29] M. Doran, S. Stern, E. Thommes, Baryon Acoustic Oscillations and Dynamical Dark Energy, JCAP 0704, 015 (2007) [arXiv:astro-ph/0609075]
  * [30] E.V. Linder, G. Robbers, JCAP 0806, 004 (2008) [arXiv:0803.2877]
  * [31] L. Knox, M. Millea, The Hubble Hunter’s Guide, Phys. Rev. D 101, 043533 (2020) [arXiv:1908.03663]
  * [32] J.C. Hill, E. McDonough, M.W. Toomey, S. Alexander, Early Dark Energy Does Not Restore Cosmological Concordance, Phys. Rev. D 102, 043507 (2020) [arXiv:2003.07355]
  * [33] K. Jedamzik, L. Pogosian, G-B. Zhao, Why reducing the cosmic sound horizon can not fully resolve the Hubble tension, arXiv:2010.04158
  * [34] V.I. Sabla, R.R. Caldwell, No $H_{0}$ Assistance from Assisted Quintessence, arXiv:2103.04999
  * [35] L. Parker, A. Raval, New quantum aspects of a vacuum-dominated universe, Phys. Rev. D 62, 083503 (2000), [Erratum: Phys. Rev. D 67, 029903 (2003)], arXiv:gr-qc/0003103
  * [36] L. Parker, D.A.T. Vanzella, Acceleration of the universe, vacuum metamorphosis, and the large-time asymptotic form of the heat kernel, Phys. Rev. D 69, 104009 (2004) [arXiv:gr-qc/0312108]
  * [37] R.R. Caldwell, W. Komp, L. Parker, D.A.T. Vanzella, A Sudden Gravitational Transition, Phys. Rev. D 73, 023513 (2006) [arXiv:astro-ph/0507622]
  * [38] E. DiValentino, E.V. Linder, A. Melchiorri, $H_{0}$ Ex Machina: Vacuum Metamorphosis and Beyond $H_{0}$, Phys. Dark Univ. 30, 100733 (2020) [arXiv:2006.16291]
  * [39] G. Benevento, W. Hu, M. Raveri, Can Late Dark Energy Transitions Raise the Hubble constant?, Phys. Rev. D 101, 103517 (2020) [arXiv:2002.11707]
  * [40] D. Camarena, V. Marra, On the use of the local prior on the absolute magnitude of Type Ia supernovae in cosmological inference, MNRAS (2021) [arXiv:2101.08641]
  * [41] G. Efstathiou, To $H_{0}$ or not to $H_{0}$?, arXiv:2103.08723
  * [42] G. Alestas, L. Perivolaropoulos, Late time approaches to the Hubble tension deforming $H(z)$, worsen the growth tension, MNRAS (2021) [arXiv:2103.04045]
  * [43] E. DiValentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D.F. Mota, A.G. Riess, J. Silk, In the Realm of the Hubble tension? a Review of Solutions, arXiv:2103.01183
  * [44] E.V. Linder, No Slip Gravity, JCAP 1803, 005 (2018) [arXiv:1801.01503]
  * [45] E.V. Linder, Limited Modified Gravity, JCAP 2010, 042 (2020) [arXiv:2003.10453]


#  A Light-Weight Multi-Objective Asynchronous Hyper-Parameter Optimizer

Gabriel Maher email<EMAIL_ADDRESS>Stephen Boyd  Mykel
Kochenderfer  Cristian Matache  Alex Ulitsky  Slava Yukhymuk  Leonid Kopman

(BlackRock AI Lab, )

###### Abstract

We describe a light-weight yet performant system for hyper-parameter
optimization that approximately minimizes an overall scalar cost function that
is obtained by combining multiple performance objectives using a target-
priority-limit scalarizer. It also supports a trade-off mode, where the goal
is to find an appropriate trade-off among objectives by interacting with the
user. We focus on the common scenario where there are on the order of tens of
hyper-parameters, each with various attributes such as a range of continuous
values, or a finite list of values, and whether it should be treated on a
linear or logarithmic scale. The system supports multiple asynchronous
simulations and is robust to simulation stragglers and failures.

##  1 Hyper-parameters, objectives, and cost

##### Hyper-parameters.

We let $x\in\mathcal{X}\subset{\mbox{\bf R}}^{n}$ denote the vector of hyper-
parameters that we will choose from the set of valid choices $\mathcal{X}$, to
be discussed in section 3). We refer to $x_{i}$ as the $i$th hyper-parameter.
While the algorithm we describe will work in principle for any value of $n$,
our focus is on the case where $n$ is modest (e.g., $<10$).

##### Objectives.

We have $k$ objectives, denoted $f_{i}:\mathcal{X}\to{\mbox{\bf R}}$, that we
use to evaluate a choice of $x$. Each of these objectives has a sense, which
means that we either want the objective to be large (i.e., to maximize it) or
small (i.e., to minimize it). We refer to these as maximization and
minimization objectives, respectively. We use $f$ to denote the vector-valued
objective, $f(x)=(f_{1}(x),\ldots,f_{k}(x))$.

##### Overall cost.

To compare different choices of hyper-parameters, we scalarize the objectives
to obtain an overall cost which we denote $F:\mathcal{X}\to{\mbox{\bf
R}}\cup\\{\infty\\}$, which we assume we wish to minimize, with $F(x)=\infty$
meaning an unacceptable value of $x$. We express this as

$F(x)=\phi(f(x))=\phi(f_{1}(x),\ldots,f_{k}(x)),$

where $\phi:{\mbox{\bf R}}^{k}\to{\mbox{\bf R}}\cup\\{\infty\\}$ is the
scalarization function, which is monotone nondecreasing in each argument
corresponding to a minimization objective, and monotone nonincreasing in each
argument corresponding to a maximization objective. (Roughly speaking, this
means that $F$ decreases when a minimization objective decreases or when a
maximization objective increases.)

A very simple scalarizer is $\phi(u)=w^{T}u$, where $w\in{\mbox{\bf R}}^{n}$
is a set of weights, with $w_{i}\geq 0$ for minimization objectives and
$w_{i}\leq 0$ for maximization objectives. In this case, $F(x)$ is a
traditional weighted sum of the objectives.

##### Target-priority-limit scalarizer.

The target-priority-limit scalarizer is a bit more complex than the weighted-
sum scalarizer, while still being expressive and interpretable. It is
separable, i.e., it has the form

$\phi(u)=\sum_{i=1}^{k}\phi_{i}(u_{i}),$

where $\phi_{i}:{\mbox{\bf R}}\to{\mbox{\bf R}}\cup\\{\infty\\}$. Each
$\phi_{i}$ is characterized by three numbers: a _target_ value $T_{i}$, a
_priority_ $P_{i}>0$, and a _limit_ $L_{i}$. For an objective that we want to
minimize, we require that $T_{i}\leq L_{i}$, and have

$\phi_{i}(u_{i})=\left\\{\begin{array}[]{ll}0&u_{i}\leq T_{i}\\\
P_{i}\frac{u_{i}-T_{i}}{L_{i}-T_{i}}&T_{i}\leq u_{i}\leq L_{i}\\\
\infty&u_{i}>L_{i}.\end{array}\right.$

We note that $\phi_{i}(u_{i})\geq 0$, taking its smallest value $0$ only when
the objective is less than or equal to the target value. When the objective
exceeds the target value but not the limit, $\phi_{i}(u_{i})$ is proportional
to the fractional distance between these two values, with slope $P_{i}$. If
$u_{i}$ exceeds the limit, we have $\phi_{i}(u_{i})=\infty$. In other words,
the target is the value of the objective below which we are equally satisfied,
and the limit specifies the maximum allowed value. The priority sets the slope
in between.

For an objective we wish to maximize, we require $L_{i}\leq T_{i}$ and reverse
the function as

$\phi_{i}(u_{i})=\left\\{\begin{array}[]{ll}0&u_{i}\geq T_{i}\\\
P_{i}\frac{T_{i}-u_{i}}{T_{i}-L_{i}}&L_{i}\leq u_{i}\leq T_{i}\\\
\infty&u_{i}<L_{i}.\end{array}\right.$

Figure 1 shows two examples of $\phi_{i}$, one for a minimization objective
and one for a maximization objective.

$-5$$0$$5$$10$$15$$20$$0$$0.5$$1$$1.5$$2$$u_{i}$$\phi_{i}(u_{i})$minimize
($T=-3,L=3,P=2$)maximize ($T=5,L=-3,P=0.5$) Figure 1: An example of the
target-priority-limit scalarizer, with one minimization objective and one
maximization objective. The overall cost is
$\phi_{1}(f_{1}(x))+\phi_{2}(f_{2}(x))$.

##  2 Hyper-parameter optimizer

##### Simulation.

We assume that the objectives $f_{1},\ldots,f_{k}$ are given as black-box
oracles that we can evaluate. We cannot evaluate gradients of $f_{i}$; in
fact, $f_{i}$ is not required to be differentiable. We refer to evaluating the
objectives $f(x)=(f_{1}(x),\ldots,f_{k}(x))$ as a _simulation_ , which can
involve substantial computation. We assume that simulations can be carried out
by a set of workers or agents.

##### Hyper-parameter optimizer.

Our goal is to find a good value of the hyper-parameters, i.e., to find
$x\in\mathcal{X}$ that approximately minimizes $F(x)$, using only simulations.
Thus we approximately solve the hyper-parameter optimization problem

$\begin{array}[]{ll}\mbox{minimize}&F(x)\\\ \mbox{subject
to}&x\in\mathcal{X},\end{array}$

with variable $x$, using only evaluations of $F$, i.e., simulations.

##### Trade-off mode.

As a variation on scalarization of the objectives, we can also search for
Pareto optimal points with respect to two or three of the objectives. Our goal
is to explore the trade-off curve (for two objectives) or trade-off surface
(three objectives). This mode is described in section 5. Until then, we
consider the scalarized hyper-parameter problem described above.

##### Hyper-parameter optimizer state.

Our hyper-parameter optimizer is a centralized entity that keeps track of
previous simulations, i.e., a list

$(x^{i},f(x^{i})),\quad i=1,\ldots,K,$

where $K$ is the number of simulations carried out so far,
$x^{i}\in\mathcal{X}$, and the superscript $i$ indexes the simulations. This
list is the state of the optimizer. For each of these simulations we can
evaluate the overall cost $F(x^{i})=\phi(f(x^{i}))$, and sort the list of
previous simulations from best (i.e., smallest $F(x^{i})$) to worst (i.e.,
largest $F(x^{i})$). The best $r$ simulations can be displayed as a leader-
board.

##### Interaction with workers.

The optimizer interacts with multiple workers (such as CPUs, processes, or
threads) that carry out simulations asynchronously. When a worker is available
to carry out a simulation, it queries the optimizer to obtain a recommended
value of $x$ to simulate. When a worker finishes a simulation, it reports
values of the $n$ hyper-parameters and the $k$ objectives found, which is
added to the list of previous simulations. We do not assume that the workers
are reliable; that is, we assume they can take a long time to evaluate the
objectives for a given value of $x$, or even fail, i.e., never report back.

##  3 Parameter attributes

Each hyper-parameter can have a number of attributes, which we describe here.
(Not all combinations of these attributes make sense.) These collectively give
the set of allowed values $\mathcal{X}$.

##### Range.

Each hyper-parameter has a minimum and maximum allowed value, denoted $l_{i}$
and $u_{i}$, with $l_{i}<u_{i}$.

##### Linear or logarithmic scale.

Each hyper-parameter has either a linear or logarithmic scale. If it is
logarithmic, we require that $l_{i}>0$. The choice will affect how values are
treated internally, and the meaning of other attributes such as gridded
values, described below.

##### Grid and finite set of values.

We can require that a hyper-parameter take on only a specified finite set of
values. For example, we can require that $x_{i}\in\\{0,1\\}$, i.e., that the
$i$th hyper-parameter is Boolean. A common special case is a grid of $N$
values, given by

$x_{i}\in\\{l_{i},l_{i}+\delta,\ldots,l_{i}+(N-1)\delta\\},\quad\delta=\frac{u_{i}-l_{i}}{N-1},$

when $x_{i}$ has linear scale, or

$x_{i}\in\\{l_{i},\gamma
l_{i},\ldots,\gamma^{N-1}l_{i}\\},\quad\gamma=(u_{i}/l_{i})^{1/(N-1)},$

when $x_{i}$ has logarithmic scale. Another common specification of a finite
set of values is that $x_{i}$ should be an integer, i.e.,

$x_{i}\in\\{\lceil l_{i}\rceil,\lceil l_{i}\rceil+1,\ldots,\lfloor
u_{i}\rfloor\\}.$

##  4 Method

###  4.1 Standarized hyper-parameters

##### Standardizing the hyper-parameters.

We first transform each $x_{i}$ to $z_{i}$, a standardized hyper-parameter
value in $[0,1]$. For a linear scale hyper-parameter, we use

$z_{i}=\frac{x_{i}-l_{i}}{u_{i}-l_{i}},$

and for a logarithmic scale hyper-parameter,

$z_{i}=\frac{\log x_{i}-\log l_{i}}{\log u_{i}-\log l_{i}}.$

The standardized hyper-parameter $z$ lies in $[0,1]^{n}$. We let $\mathcal{Z}$
correspond to $\mathcal{X}$ after standardization. Thus $z\in\mathcal{Z}$
means that the corresponding $x$ satisfies $x\in\mathcal{X}$.

##### Un-standardizing.

Our method will suggest values of $z\in{\mbox{\bf R}}^{n}$, which are
transformed back to $x$ by first projecting onto $\mathcal{Z}$, and then
inverting the standardization mappings. Given the suggestion $z\in[0,1]^{n}$,
we first form

$\tilde{z}=\Pi_{\mathcal{Z}}(z),$

where $\Pi_{\mathcal{Z}}$ is projection onto $\mathcal{Z}$. This projection
involves clipping the entries $z_{i}$ to $[0,1]$ and then rounding to the
nearest valid value of $z_{i}$. We then un-standardize using

$x_{i}=l_{i}+z_{i}(u_{i}-l_{i})$

for linear scale hyper-parameters and

$x_{i}=\exp\left(\log l_{i}+z_{i}(\log u_{i}-\log l_{i})\right)$

for logarithmic scale hyper-parameters.

###  4.2 Sampling and exploration strategy

When a worker requests a hyper-parameter value to simulate, the optimizer
generates $z$ randomly from some distribution, and then un-standardizes it as
described above to obtain a valid suggestion $x\in\mathcal{X}$. We now explain
how we represent and update this distribution.

##### Initial exploration phase.

We let $K$ denote the total number of simulations carried out so far. For $K$
less than a given threshold, we simply sample $z$ from a uniform distribution
on $[0,1]^{n}$. Our default threshold is

$\min\\{\lfloor S/5\rfloor,50+2n\\},$

where $S$ is the intended number of simulations. Once $K$ equals the
threshold, we switch to a more sophisticated distribution that depends on the
state of the optimizer, i.e., the table of previous simulations.

##### Elite points.

We work with the elite points, which is the set of the $r$ best values we have
found so far, judged in terms of overall cost. We take $r$ to be some fraction
$\eta$ (e.g., $0.2$) of $K$, the total number of simulations we have carried
out. In other words, the elite points are the ones with cost at or below the
$\eta$-quantile among the simulations.

##### Statistical model of elite points.

We fit a Gaussian mixture model of the standardized versions of the elite
points. This statistical model is updated whenever the elite set changes,
which can occur when new simulation results arrive, or when the scalarization
changes. When a worker requests a new hyper-parameter value, we sample from
this distribution, project onto $\mathcal{Z}$, and then un-standardize.

##### Low discrepancy sample points.

Instead of sampling from a uniform distribution for the initial phase, we
instead use Sobol points [1], which are more space filling than random
samples. For the Gaussian mixture model, we may also sample from points
derived from a Sobol sequence. To generate these Gauss-Sobol points, we find a
Sobol sequence in ${\mbox{\bf R}}^{2}$, and then transform these in the
standard way to obtain points in ${\mbox{\bf R}}^{2}$ that would have had an
$\mathcal{N}(0,I)$ distribution had the points been uniform. These, in turn,
are transformed by an appropriate affine mapping to obtain the desired mean
and covariance.

##  5 Trade-off mode

This tool supports multi-objective optimization in a trade-off mode. In this
case, we do not scalarize the objectives, but instead indicate two or three
for which we seek Pareto optimal points. For all other objectives we can
maintain limits, i.e., we still consider only choices for which all metrics
are within the limits, but we do not set priorities or target values. Our goal
is to explore the trade-off curve (for two objectives) or trade-off surface
(for three objectives).

To handle this trade-off exploration, we only need to change how we
characterize previous simulations as elite or not. Similar to [2] [2], we
define the level 1 simulations as those that are Pareto optimal with respect
to the given two or three objectives. Then we define the level 2 simulations
to be those that are Pareto optimal after the level 1 simulations are removed.
We continue this for other levels until all past simulations have been
assigned a level. To choose $r$ previous simulations, we first consider the
level 1 simulations. If $r$ is smaller than or equal to the number of level 1
simulations, we randomly choose $r$ of them. Otherwise we label all level 1
simulations as elite and move to level 2 simulations to find more elite
points, repeating this process. In other words, we sort the simulations
according to their Pareto level instead of a scalarized overall objective
value; otherwise, the search algorithm is the same.

##  6 The HOLA software package

We implemented the method described above in the Python package HOLA (Hyper-
parameter Optimizer: Light-weight and Asynchronous). We describe HOLA here.

##### Objective and scalarizer configuration.

The objective names, targets, limits and priorities are specified as a key-
value dictionary, for example

    objectives_config = {
        "accuracy":{
            "target": 1.0,
            "limit": 0.0,
            "priority": 2.0
        },
        "abs_error":{
            "target": 0,
            "limit": 1000,
            "priority": 0.5
        }
    }

##### Hyper-parameter search space configuration.

The search-space for each hyper-parameter is specified with a few options
supplied to our software as a simple key-value dictionary. First, users must
specify a `min` and `max` value. Thereafter, the scale can be set to either
linear or logarithmic and defaults to linear if no choice is made. Users can
also elect to sample from an equally spaced grid of points between `min` and
`max`. Note that the grid will be equally spaced on a log scale if the scale
is set to logarithmic. Finally hyper-parameter values can be sampled from a
set of discrete values as well. An example hyper-parameter search-space
configuration could look like

    params_config = {
        "n_estimators":{
            "min": 10,
            "max": 1000,
            "param_type": "int",
            "scale": "log"
            "grid": 10
        },
        "max_depth":{
            "values": [1, 3, 5, 7]
        },
        "learning_rate":{
            "min": 1e-4,
            "max": 1.0,
            "scale": "log"
        },
        "subsample":{
            "min": 0.2,
            "max": 1.0
        }
    }

##### Running hyper-parameter optimization in local mode.

HOLA exposes a `tune` function that provides a straightforward interface to
run hyper-parameter optimization sessions on a single machine (but still
utilizing all available processors). Users simply pass the objectives, the
configuration dictionaries, as well as a stopping criterion called `num_runs`
which corresponds to the number of simulations made until the best hyper-
parameters are reported. Calling

`tune(func, params_config, objectives_config, num_runs=100, n_jobs=2)`

spawns `n_jobs` processes that request hyper-parameter suggestions from HOLA
and execute `func` until `num_runs` simulations are complete. If `n_jobs=-1`,
then HOLA will detect the number of available processors and spawn that many
processes.

Users connect their simulation to HOLA using the `func` argument. The inputs
to `func` are expected to be the same as the keys in the hyper-parameter
search-space configuration dictionary. `func` is expected to return a key-
value dictionary with objective names and corresponding objective values as
keys and values respectively. In other words, `func` should accept the hyper-
parameters as input, run a simulation and return a key-value dictionary of the
resulting objective values.

The `tune` function returns an underlying `Tuner` object which is responsible
for running the optimizer and maintaining the state. While running
simulations, the `Tuner` object will simultaneously sort the leader-board of
simulation results and provide optimized hyper-parameter suggestions. Calling
`Tuner.get_best_params()` will return the current best hyper-parameter values
in the leader-board. Similarly, `Tuner.get_best_scores()` will return the best
objective and scalarized score values.

##### Saving and restoring hyper-parameter optimization sessions.

Our package stores completed simulation results in a leader-board that can be
saved to a dataframe at any time. A hyper-parameter optimization session can
be restored by loading a saved leader-board dataframe, supplying an objectives
configuration dictionary, hyper-parameter search-space configuration
dictionary, and the simulation function. It is thus straightforward to save
and restore hyper-parameter optimization sessions using our package.

##### HOLA hyper-parameter optimization server.

We provide a server implementation of the HOLA package such that the workers
can be implemented and run completely independently, for example across
multiple machines or in a different programming language. The HOLA server
exposes a simple HTTP interface through which workers can request hyper-
parameter samples and report simulation results. Note that we assume a single
server instance will be responsible for a single family of simulation results,
that is all simulation results reported by workers correspond to the same
hyper-parameter and objectives configuration. Futhermore the leaderboard of
hyper-parameters and simulation results can be viewed in a web browser in
real-time. The server is launched by specifiying a directory containing a JSON
hyper-parameter and objective configuration file named `hola_params.json` and
`hola_objectives.json` respectively. Optionally, if the directory contains a
saved leaderboard csv file named `hola_results.csv`, the leaderboard will be
restored and the server will resume from that point on. By default the server
will run on the `localhost` address on port 8675, but an alternative port can
be specified as a command-line argument.

We expose the following routes through the HOLA server for workers to request
hyper-parameter samples and report simulations results

  * •

`/`

    * –

`GET`: Will return a HTML page with the current leaderboard results, viewable
in a browser.

  * •

`/report_request`

    * –

`GET`: Will return a hyper-parameter sample as a JSON response with the keys
being the hyper-parameter names and the values being their sampled values.

    * –

`POST`

      * *

`request`: The request can be empty or optionally a JSON dictionary with keys
`params` and `objectives`. The `params` should be a dictionary of hyper-
parameter names and values. The `objectives` should be a dictionary of
objective names and values.

      * *

`response`: If the request is empty, will simply return a JSON hyper-parameter
sample dictionary. If the request is non-empty, the hyper-parameter sample and
simulation result will be recorded in the leaderboard and a JSON response with
a new hyper-parameter sample will be returned.

  * •

`/param`

    * –

`GET`: Will return a key-value dictionary of the best hyperparameters seen so
far.

  * •

`/experiment`

    * –

`GET`: Will return a JSON response dictionary with keys `params` and
`objectives` containing the hyper-parameter and objective configuration
dictionaries, respectively.

##  7 Prior work

A simple approach to hyper-parameter optimization known as full-factorial
evaluation is to discretize the various parameters values and try out all
possible combinations of discrete values, corresponding to points on a grid.
Alternatively hyper-parameter values could also be sampled from the hyper-
parameter search space uniformly at random and evaluated [3]. Because the
gradient of the simulation results with respect to the hyper-parameters is
typically unavailable, gradient-free optimization methods can also be used for
hyper-parameter optimization [4]. Search and exploration methods using
Gaussian Processes to model the hyper-parameter distribution have also been
applied [5]. Our method is a modified form of the cross-entropy method [6].
There are many open-source software packages dedicated to hyper-parameter
optimization such as Spearmint [7], Optuna [8], Tune [9], HyperOpt [10], and
SigOpt [11].

## References

  * [1] I.. Sobol  “On the Distribution of Points in a Cube and the Approximate Evaluation of Integrals”  In _USSR Computational Mathematics and Mathematical Physics_ 7.4, 1967, pp. 86–112
  * [2] K. Deb, A. Pratap, S. Agarwal and T. Meyarivan  “A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II”  In _IEEE Transactions on Evolutionary Computation_ 6.2, 2002, pp. 182–197
  * [3] James Bergstra and Yoshua Bengio  “Random Search for Hyper-Parameter Optimization” https://www.jmlr.org/papers/volume13/bergstra12a/bergstra12a In _Journal of Machine Learning Research_ , 2012
  * [4] Mykel J. Kochenderfer and Tim A. Wheeler  “Algorithms for Optimization”  MIT Press, 2019
  * [5] James Bergstra, Rémi Bardenet, Yoshua Bengio and Balázs Kégi  “Algorithms for Hyper-Parameter Optimization” https://proceedings.neurips.cc/paper/2011/file/86e8f7ab32cfd12577bc2619bc635690-Paper.pdf In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2011
  * [6] Reuven Y. Rubinstein and Dirk P. Kroese  “The Cross-Entopy Method”  Springer, 2004
  * [7] Jasper Snoek, Hugo Larochelle and Ryan P. Adams  “Practical Bayesian Optimization of Machine Learning Algorithms”  In _Advances in Neural Information Processing Systems (NIPS)_ , 2012
  * [8] Takuya Akiba et al.  “Optuna: A Next-generation Hyperparameter Optimization Framework”  In _Conference on Knowledge Discovery and Data Mining_ , 2019
  * [9] Richard Liaw et al.  “Tune: A Research Platform for Distributed Model Selection and Training”  In _arXiv preprint arXiv:1807.05118_ , 2018
  * [10] J. Bergstra, D. Yamins and D.. Cox  “Making a Science of Model Search: Hyperparameter Optimization in Hundreds of Dimensions for Vision Architectures”  In _International Conference on Machine Learning (ICML)_ , 2013
  * [11] Ian Dewancker et al.  “Evaluation System for a Bayesian Optimization Service”  In _arXiv preprint arXiv:1605.06170_ , 2016
  * [12] J.. Nelder and R. Mead  “A Simplex Method for Function Minimization”  In _The Computer Journal_ 7.4, 1965, pp. 308–313
  * [13] R. Hooke and T.. Jeeves  “Direct Search Solution of Numerical and Statistical Problems”  In _Journal of the ACM (JACM)_ 8.2, 1961, pp. 212–229
  * [14] D.. Goldberg  “Genetic Algorithms in Search”  Addison-Wesley, 1989
  * [15] J. Kennedy, R.. Eberhart and Y. Shi  “Swarm Intelligence”  Morgan Kaufmann, 2001

##  Appendix A Algorithm evaluation

This paper introduces not only a general multi-objective optimization
framework that could work with several optimization algorithms but also a new
such algorithm. To understand how our HOLA algorithm performs relative to
other well-known algorithms, we studied the quality of the best values found
on different benchmarks problems with different numbers of simulations.

###  A.1 Optimization algorithms

We compared our algorithm with a few other standard algorithms.

  1. 1.

_Random Search._ The algorithm suggests new hyper-parameters by sampling
random values for each of them.

  2. 2.

_Sobol Search._ Instead of sampling randomly, this algorithm uses a Sobol
sequence.

  3. 3.

_Random Search X2._ This algorithm is the same as Random Search but with twice
the budget (i.e., twice as many iterations).

  4. 4.

_Iterative Grid Refinement (IGR)._ This algorithm starts with a grid search
with loosely spaced lattices, finds the best parameters, and repeats the grid
search in smaller subspaces around the best parameters found so far. It
tightens the spacing between lattices in these subspaces until a stopping
criterion is met.

  5. 5.

_Nelder Mead._ A method that incrementally updates a set of points that
defines a simplex according to a set of rules [12].

  6. 6.

_Hooke Jeeves._ A method that greedily explores points some step size along
the various coordinate axes. It reduces the step size if an improvement around
a point cannot be found [13].

  7. 7.

_Tree-structured Parzen Estimators (TPE)._ A Bayesian optimization technique
[5].

  8. 8.

_HOLA._ This is our modified version of the cross-entropy method with Gaussian
mixture models.

  9. 9.

_Genetic Algorithms (GA)._ This algorithm evolves a population of points using
crossover and mutation operations [14].

  10. 10.

_Particle Swarm Optimization (PSO)._ This algorithm moves particles
corresponding to points in the design space using various forces that aim to
bring the particles close to a global optimum [15].

We used the hyperopt library for TPE, Random Search, and Random Search X2. For
Hooke Jeeves, Nelder Mead, GA, and PSO, we used the pymoo library. For HOLA,
IGR, and Sobol Search, we used our own implementation.

When running an optimizer, the user typically specifies when the algorithm
will terminate. In our experiments, we used the number of iterations as our
stopping criterion. Most of these algorithms are “informed” in the sense that
they use previous values to make better informed parameter suggestions in the
next iteration.

###  A.2 Benchmarks

Here we describe what functions we used for benchmarking our optimization
algorithms. We ranged the dimensionality of the input domains of these
functions between 1D and 7D .

In terms of real-world scenarios, we maximized the R-squared (by minimizing
its opposite) of GradientBoostingRegressor from scikit-learn on the built-in
scikit-learn Diabetes Regression dataset with 4 hyper-parameters: number of
estimators, maximum depth, learning rate and subsample.

Moreover, we experimented with several closed-form hard-to-optimize functions:
Ackley (2D, 5D, 7D), Branin (2D), Bukin 6 (2D), Cross In Tray (2D), Drop Wave
(2D), Egg Holder (2D), Forrester (1D), Holder Table (2D), Levy 13 (2D),
Rastrigin (2D, 5D, 7D), Schwefel (2D, 5D, 7D) and Six Hump Camel (2D) .

The difficulty in optimizing a function is determined by its modality
(unimodal, bimodal, multimodal), smoothness and dimensionality. For the
functions that support an arbitrary number of dimensions (e.g., Rastrigin),
the higher the dimensionality the harder the problem is to optimize. Also, the
smoother the function, the easier it should be for an informed optimizer to
optimize it. Most of the functions we chose are multimodal, and some are maybe
quite exaggerated for what would occur in real-world, optimization problem,
but we wanted to see how far HOLA can go by presenting it challenging
problems. To balance this, some chosen functions like Branin, Forrester, or
Six Hump Camel have fewer local minima.

###  A.3 Experiments and Results

We ran many times each of the enumerated algorithms on the mentioned
benchmarks ranging the number of iterations.

####  A.3.1 Real-world example

For the gradient boosted regressors, we ran at least 50 times each algorithm
at 25, 50, 75, 100 and 200 iterations. Figure 2 shows the distribution of the
best values found by each optimizer at the different numbers of iterations in
a box plot. The median is shown as a solid line while and the mean as a dotted
one. Since this was implemented as a minimization problem, the lower is
better. The HOLA algorithm shows that despite its simplicity it is very
competitive when compared to the other algorithms. Note, that it performs
significantly better than TPE which is a method that belongs to the same
Bayesian class.

Figure 2: Best values found for Gradient Boosted Regressor on the Diabetes
dataset

####  A.3.2 Generalization on Synthetic Benchmarks

With this encouraging real-world results, we ran further experiments to
understand how well HOLA generalizes. Therefore, we turned to the synthetic
hard-to-optimize functions described in section A.2.

Since the performance of algorithms is directly proportional with the number
of iterations, for the smaller numbers of iterations (25, 50 and 75), we used
”easier” benchmarks: maximum 2 dimensional. That is: Ackley (2D), Branin (2D),
Bukin 6 (2D), Cross In Tray (2D), Drop Wave (2D), Egg Holder (2D), Forrester
(1D), Holder Table (2D), Levy 13 (2D), Rastrigin (2D), Schwefel (2D), Six Hump
Camel (2D) . For higher number of iterations, we used ”harder” benchmarks:
also included 5 and 5 dimensional functions and eliminated the “easier”
functions where all optimizers were optimizing nearly equally well: Forrester,
Branin, Six Hump Camel, Cross In Tray. So, for the higher number of iterations
we are left with Ackley (2D, 5D, 7D), Bukin 6 (2D), Drop Wave (2D), Egg Holder
(2D), Holder Table (2D), Levy 13 (2D), Rastrigin (2D, 5D, 7D), Schwefel (2D,
5D, 7D).

Again we ran each optimizer on each benchmark for each number of iterations at
least 50 times. In total, we collected ran 10 optimization algorithms which
tried 29668300 combinations of hyper-parameters on 18 multimodal functions
(i.e., benchmarks).

To compare across benchmarks we min-max normalized these values and computed,
first the mean per number of iterations, optimizer and benchmark, and then we
took again the mean of these means per number of iterations and optimizer.
That is, we can observe the performance of the optimizers for each number of
iterations, across all benchmarks by plotting this ”mean of means of
normalized best values found by simulations”.

Figure 3: Average performance per optimizer across hard-to-optimize functions
Figure 4: Average performance per optimizer across hard-to-optimize functions

In figures 3 and 4, we show these means. Since all problems are minizimations
lower is better. These plots reinforce the existing literature as:

  * •

The more iterations the better an algorithm performs.

  * •

Genetic Algorithms and PSO perform better at a high number of iterations.

  * •

Random Search X2 is better than Sobol Search which in turn is better than
Random Search.

Based on this further analysis, we can conclude that HOLA is performant in a
more general context of multimodal hard-to-optimize functions. Random Search
X2 is so cheap in terms of resources and simple to implement, but HOLA
significantly outperforms it with half the resources. Moreover, HOLA dominates
the popular TPE algorithm (hyperopt), the latter only catching up with HOLA at
1000 iterations. When compared to the other algorithms, HOLA is quite simple
but nonetheless powerful. Similarly, the Iterative Grid Refinement (IGR)
algorithm is also very simple and quite performant.