#  On convergence of volume of level sets of stationary smooth Gaussian fields

Dmitry Beliaev 111Mathematical Institute, University of Oxford, UK  Akshay
Hegde 11footnotemark: 1

###### Abstract

We prove convergence of Hausdorff measure of level sets of smooth Gaussian
fields when the levels converge. Given two coupled stationary fields
$f_{1},f_{2}$ , we estimate the difference of Hausdorff measure of level sets
in expectation, in terms of $C^{2}$-fluctuations of the field $F=f_{1}-f_{2}$.
The main idea in the proof is to represent difference in volume as an integral
of mean curvature using the divergence theorem. This approach is different
from using Kac-Rice type formula as main tool in the analysis.

††Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>††2020 Mathematics
Subject Classification. 60G60, 60G15, 53A07. ††Keywords and phrases. Gaussian
fields, level sets.

##  1 Introduction

Smooth Gaussian fields appear naturally in some areas of mathematics and
physics. Random plane waves are conjectured to describe local behaviour of
particles in chaotic domains, to give an example [5]. It also interacts well
with other topics such as percolation theory [3]. Random fields have found
some applications in areas as diverse as oceanography [7], cosmology [2],
medical imaging [11].

Studying geometrical and topological properties of the field, especially of
level/excursion sets of the field is of great interest. Particularly,
functionals such as volume of level sets, number connected components of level
sets are well studied (see [10],[9]). In problems involving Gaussian fields,
sometimes one needs to compare two fields, say by coupling them, when their
laws are close. Comparing geometric observables are of particular interest. We
show that, with probability close to one, difference in Hausdorff measures of
_nodal sets_ (i.e. the zero sets) of coupled fields with ‘close’ laws is
small. The main idea in the proof is to represent difference in volumes of
level sets as an integral of mean curvature of the hypersurface using the
divergence theorem. This representation is classical in Riemannian geometry
and has been used extensively in study of minimal surfaces [6, Chapter 1]. The
novelty is to get an average estimate of the difference in volumes in the
context of Gaussian fields. Also, we don’t rely on Kac-Rice (or any other
variation of co-area formula) for the analysis of volume of level sets, which
is a standard tool in this topic. As a by product, we give an explicit formula
for the mean curvature of level sets at a given level. We believe that proving
convergence in distribution of Hausdorff measure of level sets can be done by
following the proof idea of Kac-Rice as presented in [1, Theorem 6.2]. But it
might not be as straight forward as our proof, and proving other convergences
might require some new ideas.

##  2 Results

In this article, we consider smooth Gaussian fields
$f:\mathbb{R}^{d}\to\mathbb{R}$ with mild non-degeneracy conditions, of fixed
dimension $d\geq 2$. Call a field _stationary_ if the covariance kernel
$K(x,y)=\mathbb{E}[f(x)f(y)]$ is translation invariant. Now for stationary
fields, the kernel $K$ is a Fourier transform of a finite symmetric Borel
measure $\rho$, called spectral measure (normalised to probability measure in
this article).

Fix a domain $D=[-R,R]^{d}\subset\mathbb{R}^{d}$. Consider two $C^{2}$-smooth
Gaussian fields $f_{1},f_{2}:\mathbb{R}^{d}\to\mathbb{R}$ and a coupling of
the fields $f_{1},f_{2}$, by abuse of notation, such that $F=f_{1}-f_{2}$ has
the $C^{2}$-fluctuations

$\sigma^{2}_{D}:=\sup_{x\in D}\sup_{|\alpha|\leq
2}\text{Var}[\partial^{\alpha}F(x)].$

######  Assumptions 2.1.

Assume that the fields $f_{1},f_{2}$ are

  1. 1.

stationary, $C^{2}$-smooth a.s.

  2. 2.

non-degenerate, i.e. $(f_{i},\nabla f_{i})$ has density in $\mathbb{R}^{d+1}$
for $i=1,2$.

  3. 3.

Morse functions a.s.

Let $\mathcal{L}^{n}$ denote $n$-dimensional Lebesgue measure. Let
$\mathcal{H}^{n}$ denote the $n$\- dimensional Hausdorff measure, which is
scaled so that $\mathcal{H}^{n}([0,1]^{n})=\mathcal{L}^{n}([0,1]^{n})$. Note
that by Bulinskaya lemma (see [8, section 5.3]), a.s. nodal sets are sub-
manifolds of $\mathbb{R}^{d}$ of co-dimension one. So we interchangeably use
the terms volume and Hausdorff measure.

######  Theorem 2.2.

Let $\mathcal{H}^{d-1}(f_{i}^{-1}(a))$ denote the volume of level sets in the
domain $D$. With the setup as above, we have

$\mathbb{E}|\mathcal{H}^{d-1}(f_{1}^{-1}(0))-\mathcal{H}^{d-1}(f_{2}^{-1}(0))|\leq
C(f_{1},f_{2})(\mathcal{L}^{d}(D)\sqrt{\log R})\sigma_{D}^{1/7}$

assuming $\sigma_{D}$ is small enough (say, $\sigma_{D}<1$). Here, the
constant $C(f_{1},f_{2})$ depends only on the laws of the fields and not the
coupling.

The factor $\sqrt{\log R}$ appearing in the above theorem is from the
quantitative version of Kolmogorov’s existence theorem for smooth Gaussian
fields as stated in [8, Appendix A]. Also, the exponent $1/7$ in
$\sigma_{D}^{1/7}$ is not optimal, and can be made close to $1/4$. We believe
optimal exponent of $\sigma_{D}$ is $1$ due to cancellations in the integral
of mean curvature in the bulk.

We make some comments on the assumptions on the fields. We would like to point
out that the proof of Theorem 2.2 works for non-stationary, non-degenerate
fields with positive lower bounds on fluctuations of the field and its
derivatives, with suitable modifications. We have not used stationarity in any
serious way, but the estimations of the bounds become much easier if we assume
stationarity. Also, assumption that the fields are a.s. Morse functions is not
very restrictive and many interesting non-degenerate fields we know are Morse
functions a.s. It can be shown that stationary fields with spectral measures
containing an open sets are Morse a.s. If the field is isotropic, then also we
can show that the field is Morse a.s. In particular, random plane wave (RPW)
model and Bargmann-Fock field (on $\mathbb{R}^{d}$) are Morse a.s.

One such coupling of fields is available using coupling of white noises (see
[4]). The coupling as in [4] gives the following estimate for the fluctuations
of the field $F=f_{1}-f_{2}$. We have,

$\sigma_{D}^{2}\leq
C(R^{d}+1)\inf_{\rho\in\mathcal{P}(\rho_{1},\rho_{2})}\int(|s|^{2}+|t|^{2}+1)^{2+1}|s-t|^{2}d\rho(s,t)$

where $\mathcal{P}(\rho_{1},\rho_{2})$ is the space of all symmetric couplings
of $\rho_{1}$ and $\rho_{2}$ and $C$ is a absolute constant.

Now by the coupling techniques mentioned above, $\sigma_{D}$ can be controlled
by the transport distances between the measures in the domain (in the general
case) or by norm of differences in spectral densities (in special cases).
Let’s give an example where this is useful. Recall that the spectral measure
of random planar waves is the uniform measure on the unit circle in
$\mathbb{R}^{2}$. We can approximate this measure, in the transport distance
mentioned above, by a measure supported on finite points. This field
corresponds to a finite interference of pure sine waves. So we can
quantitative bounds on the difference of lengths.

To prove Theorem 2.2, first we study convergence of volume of level sets using
the divergence theorem. Although expressing change in volume of a hypersurface
in normal direction in terms of mean curvature is classical as previously
mentioned, we need the version as in Proposition 2.3.

######  Proposition 2.3.

Let $f:\mathbb{R}^{d}\to\mathbb{R}$ be a non-degenerate, $C^{2}$-smooth
Gaussian field which is Morse function a.s. Let $\mathcal{H}^{d-1}(f^{-1}(a))$
denote the volume of level set $f^{-1}(a)$ in $D$. Then, almost surely, we
have

$\mathcal{H}^{d-1}(f^{-1}(b))-\mathcal{H}^{d-1}(f^{-1}(a))=\iint_{D}\kappa\mathbbm{1}_{f\in[a,b]}d\text{vol}-\oint_{\partial
D}\left\langle\dfrac{\nabla f}{|\nabla
f|},\hat{\eta}\right\rangle\mathbbm{1}_{f\in[a,b]}d\text{S}$ (1)

where

$\kappa=\text{div}\left(\dfrac{\nabla f}{|\nabla f|}\right)$

is ($(d-1)$-times) the mean curvature of level set of $f$ at $x$ and
$\hat{\eta}$ is the outward unit normal to $\partial D$ . We also have

$\mathcal{H}^{d-1}(f^{-1}(b))\to\mathcal{H}^{d-1}(f^{-1}(a)),\text{ as }b\to
a$

almost surely and in $L^{1}$.

As a corollary, we get the following formula for the mean curvature of level
sets at a given level. Usually, it is hard to get such explicit formula for
general fields.

######  Corollary 2.4.

With assumptions as in Theorem 2.2 , we have

$\mathbb{E}[\kappa|f=a]=-a\mathbb{E}[|\nabla f|].$

##  3 Proofs

###### Proof of Proposition 2.3.

Note that $f$ has only finitely many critical points in $D$ a.s. We prove in
subsection 3.1 that $\kappa$ as a function on $D$ is integrable almost surely.
We also can assume that $f$ has no critical points on $\partial D$. This is
because of Bulinskaya lemma, since $\partial D$ is $(d-1)$-dimensional and for
non-degenerate, smooth Gaussian $f$ the gradient $\nabla f$ has (Gaussian)
density on $\mathbb{R}^{d}$.

Case 1: $a,b$ are regular values of $f$.

Let $R^{\prime}=D\cap f^{-1}[a,b]$ and the unit outward normal
$\hat{\eta}=-{\nabla f}/{|\nabla f|}$ on $f^{-1}(a)$, $\hat{\eta}={\nabla
f}/{|\nabla f|}$ on $f^{-1}(b)$ (assuming $a<b$), outward normal on parts of
$\partial D\cap f^{-1}(a,b)$. Assume that $R^{\prime}$ has no critical points
of $f$ and we know that $\kappa$ is continuous except at critical points of
$f$. Apply Greens formula for the function ${\nabla f}/{|\nabla f|}$ on
$R^{\prime}$, we get

$\begin{split}\int_{f^{-1}(b)\cap D}\left\langle\dfrac{\nabla f}{|\nabla
f|},\hat{\eta}\right\rangle dS+&\int_{f^{-1}(a)\cap
D}\left\langle\dfrac{\nabla f}{|\nabla f|},\hat{\eta}\right\rangle dS+\\\
&\oint_{\partial D}\left\langle\dfrac{\nabla f}{|\nabla
f|},\hat{\eta}\right\rangle\mathbbm{1}_{f\in[a,b]}dS=\iint_{R^{\prime}}\text{div}\left(\dfrac{\nabla
f}{|\nabla f|}\right)d\text{vol}.\end{split}$ (2)

But first two terms of LHS of above equation are
$\mathcal{H}^{d-1}(f^{-1}(b)),-\mathcal{H}^{d-1}(f^{-1}(a))$ respectively.
Hence we get equation (1) in this case.

If $R^{\prime}$ has critical points of $f$, then the number of critical points
has to be finite. Let $\\{x_{1},x_{2},\ldots x_{k}\\}$ be the critical points
in $R^{\prime}$. Now apply the divergence theorem to the field $\nabla
f/|\nabla f|$ on $R^{\prime}\setminus\cup_{j}B_{\delta}(x_{j})$. Letting
$\delta\to 0$, and using integrability $\kappa$ on $D$ (see subsection 3.1),
we again get equation (1).

Case 2: $a$ or $b$ (or both) are critical values of $f$.

First, let us show continuity of volume of level sets at all levels, including
at critical values of $f$. Fix a critical value $a$ of $f$. By Morse lemma,
$f$ can be made a quadratic function at a critical point by re-
parametrisation. Let $p$ be a critical point, then there is a neighborhood $U$
of $p$ and a smooth chart $(y_{1},y_{2},\ldots,y_{d})$ such that $y_{i}(p)=0$
and

$f(y)=f(p)\pm y_{1}^{2}\pm y_{2}^{2}\cdots\pm y_{d}^{2}.$

We know that the volume of level sets of quadratic functions are continuous.
So, given a critical point $p$ of $f$ at level $a$, volume of level sets of
$f$ in a neighborhood $U$ of $p$ converge when the levels converge to $a$.
When $x_{0}\in f^{-1}(a)$ is a regular point, then there exists a neighborhood
$U_{x_{0}}$ such that the volume of level sets are continuous. This follows
from the implicit function theorem. Now, using compactness of $f^{-1}(a)\cap
D$, we get that volume of level sets is continuous at any arbitrary level.

Since the number of critical values of $f$ is finite in $D$, any critical
level in $D$ can be approximated by regular levels of $f$ in $D$. Let
$\epsilon_{n}$ be a sequence converging to zero such that
$(b-\epsilon_{n}),(a+\epsilon_{n})$ are sequences of regular values of $f$. By
continuity of the volume of level sets, we have

$\mathcal{H}^{d-1}(f^{-1}(b))-\mathcal{H}^{d-1}(f^{-1}(a))=\lim_{n\to\infty}[L(b-\epsilon_{n})-L(a+\epsilon_{n})].$

Using case 1, we have the integral formula for difference of volume of level
sets. Note that

$\left\langle\dfrac{\nabla f}{|\nabla
f|},\hat{\eta}\right\rangle\mathbbm{1}_{f\in[a+\epsilon_{n},b-\epsilon_{n}]}\to\left\langle\dfrac{\nabla
f}{|\nabla f|},\hat{\eta}\right\rangle\mathbbm{1}_{f\in[a,b]},$
$\kappa\mathbbm{1}_{f\in[a+\epsilon_{n},b-\epsilon_{n}]}\to\kappa\mathbbm{1}_{f\in[a,b]}$

pointwise. Hence by the dominated convergence theorem, we have equation (1)
for case 2 as well.

We have that $\mathcal{H}^{d-1}(f^{-1}(b))\to\mathcal{H}^{d-1}(f^{-1}(a))$ as
$b\to a$ a.s. by above discussion of continuity of length w.r.t levels. We
also have
$\mathbb{E}[\mathcal{H}^{d-1}(f^{-1}(b))]\to\mathbb{E}[\mathcal{H}^{d-1}(f^{-1}(a))]$
when $b\to a$ by Kac-Rice formula. Hence, by Scheffe’s lemma, we have $L^{1}$
convergence.

∎

###### Proof of Corollary 2.4.

Take expectation to both sides of the equation (1). Switching integration and
expectation because of Fubini’s theorem, we get

$\mathbb{E}[\mathcal{H}^{d-1}(f^{-1}(b))]-\mathbb{E}[\mathcal{H}^{d-1}(f^{-1}(a))]+\oint_{\partial
D}\mathbb{E}\left[\left\langle\dfrac{\nabla f}{|\nabla
f|},\hat{\eta}\right\rangle\mathbbm{1}_{f\in[a,b]}\right]dS=\iint_{D}\mathbb{E}[\kappa\mathbbm{1}_{f\in[a,b]}]d\text{vol}.$

Now, dividing above equation by $b-a$ and taking the limit $b\to a$, we get

$-ap(a)\mathcal{L}^{d}(D)\mathbb{E}[|\nabla f|]+\oint_{\partial
D}\mathbb{E}\left[\left\langle\dfrac{\nabla f}{|\nabla
f|},\hat{\eta}\right\rangle|f=a\right]p(a)d\text{S}=\iint_{D}\mathbb{E}[\kappa|f=a]p(a)d\text{vol}$
(3)

where $p$ is the pdf of standard Gaussian random variable.

First, from stationary Kac-Rice formula, we have

$\lim_{b\to
a}\dfrac{\mathbb{E}[\mathcal{H}^{d-1}(f^{-1}(b))]-\mathbb{E}[\mathcal{H}^{d-1}(f^{-1}(a))]}{b-a}=-ap(a)\mathcal{L}^{d}(D)\mathbb{E}[|\nabla
f|].$

Next, from the continuity of the Gaussian regression formula, we get the
respective conditional expectations (see [1, Theorem 3.2] for an explanation).
Consider the expression $\mathbb{E}[\kappa\mathbbm{1}_{f\in[a,b]}]$ and write
it in the following form,

$\mathbb{E}[\kappa\mathbbm{1}_{f\in[a,b]}]=\int_{a}^{b}\mathbb{E}[\kappa|f=u]p(u)du$

Now note that $\mathbb{E}[\kappa|f=u]$ is continuous in $u$, hence
$(b-a)^{-1}\mathbb{E}[\kappa\mathbbm{1}_{f\in[a,b]}]\to\mathbb{E}[\kappa|f=a]$
as $b\to a$. By the dominated convergence theorem, we have

$\lim_{b\to
a}\dfrac{1}{b-a}\iint_{D}\mathbb{E}[\kappa\mathbbm{1}_{f\in[a,b]}]d\text{vol}=\iint_{D}\mathbb{E}[\kappa|f=a]p(a)d\text{vol}.$

A similar argument works for the claim

$\lim_{b\to a}\dfrac{1}{b-a}\oint_{\partial
D}\mathbb{E}\left[\left\langle\dfrac{\nabla f}{|\nabla
f|},\hat{\eta}\right\rangle\mathbbm{1}_{f\in[a,b]}\right]d\text{S}=\oint_{\partial
D}\mathbb{E}\left[\left\langle\dfrac{\nabla f}{|\nabla
f|},\hat{\eta}\right\rangle|f=a\right]p(a)d\text{S}.$

So, we have the equation (3).

Now, we claim that

$\oint_{\partial D}\mathbb{E}\left[\left\langle\dfrac{\nabla f}{|\nabla
f|},\hat{\eta}\right\rangle|f=a\right]d\text{S}=0.$

Since $\nabla f$ and $f$ are pointwise independent r.v. and by stationary,
integral on a $(d-1)$-dimensional slab in $\partial D$ cancels that from the
opposite slab. Again by stationarity of $\kappa$, the equation (3) reduces to
$\mathbb{E}[\kappa|f=a]p(a)=-a\mathbb{E}[|\nabla f|]p(a)$. Hence we have the
Corollary 2.4.

∎

###### Proof of Theorem 2.2.

First, observe that $\mathcal{H}^{d-1}(f^{-1}(a))\to 0$ almost surely as
$a\to\infty$ or as $a\to-\infty$, since probability that $f$ is unbounded on
$D$ is zero. Now, taking difference of equation (1) applied to $f_{1},f_{2}$
and taking $b=0$, $a\to-\infty$ we have,

$\begin{split}\mathcal{H}^{d-1}(f_{1}^{-1}(0))-\mathcal{H}^{d-1}(f_{2}^{-1}(0))=&\iint_{D}\left[\kappa_{1}\mathbbm{1}_{f_{1}\leq
0}-\kappa_{2}\mathbbm{1}_{f_{2}\leq 0}\right]d\text{vol}\\\ &-\int_{\partial
D}\left[\left\langle\dfrac{\nabla f_{1}}{|\nabla
f_{1}|},\hat{\eta}\right\rangle\mathbbm{1}_{f_{1}\leq
0}-\left\langle\dfrac{\nabla f_{2}}{|\nabla
f_{2}|},\hat{\eta}\right\rangle\mathbbm{1}_{f_{2}\leq 0}\right]dS.\end{split}$
(4)

We bound the bulk term and the boundary term of equation (4) separately.

Bulk term: First we have,

$\begin{split}\left|\int_{D}\left[\kappa_{1}\mathbbm{1}_{f_{1}\leq
0}-\kappa_{2}\mathbbm{1}_{f_{2}\leq
0}\right]d\text{vol}\right|\leq&\left|\int_{D}(\kappa_{1}-\kappa_{2})\mathbbm{1}[f_{1},f_{2}<0]d\text{vol}\right|\\\
&+\left|\int_{D}\kappa_{1}\mathbbm{1}[f_{1}f_{2}<0]d\text{vol}\right|+\left|\int_{D}\kappa_{2}\mathbbm{1}[f_{1}f_{2}<0]d\text{vol}\right|.\end{split}$
(5)

For the second term of equation (5) we show that, with probability close to
one, $\mathcal{L}^{d}(f_{1}f_{2}<0)$ is small and that integral of curvature
is bounded with high probability.

Note that $\mathbb{E}[|\kappa_{1}|^{1+\alpha}]<\infty$ for all $0<\alpha<1$
(see section 3.1). Take $\alpha=1/2$ when applying Hölder inequality in the
following computation. Given a point $x\in D$, recall that $\kappa_{1}(x)$ is
the curvature of the level set $f^{-1}(c)$, where $x\in f^{-1}(c)$, at $x$.

$\begin{split}\left|\mathbb{E}\int_{D}\kappa_{1}\mathbbm{1}[f_{1}f_{2}<0]d\text{vol}\right|&\leq\mathbb{E}\left|\int_{D}\kappa_{1}\mathbbm{1}[f_{1}f_{2}<0]d\text{vol}\right|\\\
&\leq\int_{D}\mathbb{E}|\kappa_{1}\mathbbm{1}[f_{1}f_{2}<0]|d\text{vol}\\\
&\leq(\mathbb{E}|\kappa_{1}|^{3/2})^{2/3}\int_{D}\mathbb{P}[f_{1}(x)f_{2}(x)<0]^{1/3}d\text{vol}\\\
&\leq
C_{1}\cdot\mathcal{L}^{d}(D)\sup_{D}[(\arccos(\rho(x)))^{1/3}]\end{split}$ (6)

where $\rho(x)$ is the correlation between $f_{1}(x)$ and $f_{2}(x)$,and the
constant $C_{1}$ depends only on the law of the fields. We have used the
identity (16) in the last inequality. Note that
$\arccos(x)=c_{1}\sqrt{(1-x)}+O((1-x)^{3/2})$ near $x=1$, where $c_{1}$ is a
universal constant. We have that $|1-\rho(x)|\leq\sigma^{2}_{D}/2$ for all
$x\in D$. Hence we have,

$\mathbb{E}\left[\left|\int_{D}\kappa_{1}\mathbbm{1}[f_{1}f_{2}<0]d\text{vol}\right|\right]\leq
C_{2}\mathcal{L}^{d}(D)\sigma_{D}^{1/3}$ (7)

where the constant $C_{2}$ only depends on the spectral measure.

Next, we’ll bound the term

$\mathbb{E}\left[\left|\int_{D}(\kappa_{1}-\kappa_{2})\mathbbm{1}[f_{1},f_{2}<0]d\text{vol}\right|\right].$

We split the computation into two cases: $||\nabla f_{i}||<\delta$ for one of
the $i=1,2$ and $||\nabla f_{i}||>\delta$ for both $i$’s (for some fixed
$\delta>0$).

Now,

$\begin{split}\mathbb{E}\left[\left|\int_{D}(\kappa_{1}-\kappa_{2})\mathbbm{1}[||\nabla
f_{1}||<\delta]d\text{vol}\right|\right]&\leq\int_{D}\mathbb{E}\left[|\kappa_{1}-\kappa_{2}|\mathbbm{1}[||\nabla
f_{1}||<\delta]\right]d\text{vol}\\\
&\leq(\mathbb{E}|\kappa_{1}-\kappa_{2}|^{4/3})^{3/4}\int_{D}\mathbb{P}(||\nabla
f_{1}||^{2}<\delta^{2})^{1/4}d\text{vol}\\\ &\leq
C_{3}\mathcal{L}^{d}(D)\sqrt{\delta}.\end{split}$ (8)

In the second inequality, we used the fact that curvature has $1+\alpha$
moments for $\alpha\in([0,1))$ and applied Hölder’s inequality. Observe that
$||\nabla f_{1}||^{2}$ has bounded pdf around zero, so $\mathbb{P}(||\nabla
f_{1}||^{2}<\delta^{2})=O(\delta^{2})$.

Define

$\beta:=||f_{1}-f_{2}||_{C^{2}(D)}.$

We exploit explicit representation of the curvature (14) in terms of
derivatives of the field. Given that $||\nabla f_{1}||,||\nabla
f_{2}||>\delta$ we have,

$|\kappa_{1}-\kappa_{2}|\leq\dfrac{1}{\delta^{3}}(\beta
p_{1}+\beta^{2}p_{2}+\beta^{3}p_{3})$

where $p_{i}$’s are polynomials in the first two derivatives of $f_{1}$ of
degree at most 2. Hence,

$\begin{split}\mathbb{E}\left|\int_{D}(\kappa_{1}-\kappa_{2})\mathbbm{1}[||\nabla
f_{1}||,||\nabla
f_{2}||>\delta]d\text{vol}\right|&\leq\delta^{-3}\int_{D}\mathbb{E}[(\beta
p_{1}+\beta^{2}p_{2}+\beta^{3}p_{3})]d\text{vol}.\end{split}$ (9)

Using Cauchy-Schwartz inequality and the fact that laws of the polynomials
$p_{i}$s are translation invariant, we have the following estimate,

$\mathbb{E}\left|\int_{D}(\kappa_{1}-\kappa_{2})\mathbbm{1}[||\nabla
f_{1}||,||\nabla
f_{2}||>\delta]d\text{vol}\right|\leq\dfrac{C_{4}\mathcal{L}^{d}(D)}{\delta^{3}}(\sqrt{\mathbb{E}\beta^{2}}+\sqrt{\mathbb{E}\beta^{4}}+\sqrt{\mathbb{E}\beta^{6}}).$

But by the moment estimates of $\beta$ given in [8, A.11.1], the
$\mathbb{E}\beta^{2}$ term dominates when the coupling of the fields
$f_{1},f_{2}$ close. So we have,

$\mathbb{E}\left|\int_{D}(\kappa_{1}-\kappa_{2})\mathbbm{1}[||\nabla
f_{1}||,||\nabla
f_{2}||>\delta]d\text{vol}\right|\leq\dfrac{C_{5}\mathcal{L}^{d}(D)}{\delta^{3}}(\sqrt{\mathbb{E}\beta^{2}}).$
(10)

Boundary term: We come to the boundary term of equation (4).

$\begin{split}\int_{\partial D}\left[\left\langle\dfrac{\nabla f_{1}}{|\nabla
f_{1}|},\hat{\eta}\right\rangle\mathbbm{1}_{f_{1}\leq
0}-\left\langle\dfrac{\nabla f_{2}}{|\nabla
f_{2}|},\hat{\eta}\right\rangle\mathbbm{1}_{f_{2}\leq
0}\right]dS=\int_{\partial D}\left[\left\langle\dfrac{\nabla f_{1}}{|\nabla
f_{1}|}-\dfrac{\nabla f_{2}}{|\nabla
f_{2}|},\hat{\eta}\right\rangle\mathbbm{1}[f_{1},f_{2}<0]\right]dS\\\
+\int_{\partial D}\left[\left\langle\dfrac{\nabla f_{1}}{|\nabla
f_{1}|},\hat{\eta}\right\rangle\mathbbm{1}[f_{1}<0,f_{2}>0]\right]dS+\int_{\partial
D}\left[\left\langle\dfrac{\nabla f_{2}}{|\nabla
f_{2}|},\hat{\eta}\right\rangle\mathbbm{1}[f_{2}<0,f_{1}>0]\right]dS\end{split}$
(11)

The analysis of bounds of first term of RHS of equation (11) is similar to
that of equation (8). We get that,

$\left|\mathbb{E}\int_{\partial D}\left[\left\langle\dfrac{\nabla
f_{1}}{|\nabla f_{1}|}-\dfrac{\nabla f_{2}}{|\nabla
f_{2}|},\hat{\eta}\right\rangle\mathbbm{1}[f_{1},f_{2}<0]\right]dS\right|\leq
C_{6}\mathcal{L}^{d-1}(\partial D)(\delta_{1}^{2}+\mathbb{E}\beta/\delta_{1})$
(12)

for $\delta_{1}>0$.

Now, second term of RHS is bounded by $C\cdot\mathcal{L}^{d-1}(\partial
D\cap\\{f_{1}f_{2}<0\\})$ since ${\nabla f_{1}}/{|\nabla f_{1}|}$ is unit
vector. By similar argument which lead to equation (7), we have

$\mathbb{E}\mathcal{L}^{d-1}(\partial D\cap\\{f_{1}f_{2}<0\\})\leq
C_{8}\mathcal{L}^{d-1}(\partial D)\sigma_{D}.$ (13)

This is again dominated by the quantity of RHS of equation (7).

Analysis of the final bound: We combine the bounds from (7), (8),(10), and
(12). Finally, we get

$\mathbb{E}|\mathcal{H}^{d-1}(f_{1}^{-1}(0))-\mathcal{H}^{d-1}(f_{2}^{-1}(0))|\leq
C\mathcal{L}^{d}(D)\left(\sigma_{D}^{1/3}+\sqrt{\delta}+\dfrac{\sqrt{\mathbb{E}\beta^{2}}}{\delta^{3}}+\dfrac{\delta_{1}^{2}}{R}+\dfrac{\mathbb{E}\beta}{\delta_{1}R}\right).$

Estimates from [8, A.9, A.11.1] gives us $\mathbb{E}\beta\leq
C_{1}(R)\sigma_{D}$ and $\sqrt{\mathbb{E}\beta^{2}}\leq C_{2}(R)\sigma_{D}$,
where we can show that $C_{1}(R),C_{2}(R)$ behave like $\sqrt{\log R}$.
Choosing $\delta=\sigma_{D}^{2/7},\delta_{1}=\sigma_{D}^{1/2}$, and assuming
$\sigma_{D}$ is small enough we have,

$\mathbb{E}|\mathcal{H}^{d-1}(f_{1}^{-1}(0))-\mathcal{H}^{d-1}(f_{2}^{-1}(0))|\leq
C(R)\mathcal{L}^{d}(D)\sigma_{D}^{1/7}.$

∎

###  3.1 Technical bits

Moments of curvature r.v.: We show that the $(1+\alpha)$-moments are finite,
where $0\leq\alpha<1$ for the r.v. $\kappa$ of a $C^{2}$-smooth, non-
degenerate, stationary field $f$. Observe that

$\kappa=\dfrac{|\nabla f|^{2}\text{Tr}(H(f))-\nabla fH(f)\nabla
f^{\text{T}}}{|\nabla f|^{3}}$ (14)

where $H(f)$ is the Hessian of the function $f$, by a simple algebraic
computation.

First let us prove that $\mathbb{E}[|\kappa|^{1+\alpha}]<\infty$ for $d=2$
case. The general case follows from similar computation. Observe that
$\mathbf{x}=(x_{1},x_{2},x_{3},x_{4},x_{5})=(\partial_{x}f,\partial_{y}f,\partial_{xx}f,\partial_{xy}f,\partial_{yy}f)$
is a Gaussian vector and that $(\partial_{x}f,\partial_{y}f)\text{ and
}(\partial_{xx}f,\partial_{xy}f,\partial_{yy}f)$ are independent, by
stationarity of the field $f$. Let $\Sigma$ be the covariance matrix of the
Gaussian vector $(\partial_{x}f,\partial_{y}f)$ and $\mathbb{P}_{1}$ be the
law of $(\partial_{xx}f,\partial_{xy}f,\partial_{yy}f)$. Let
$\mathbf{x}=(x_{1},x_{2})$ and $\mathbf{x}^{\prime}=(x_{3},x_{4},x_{5})$.

So,

$\begin{split}\mathbb{E}[|\kappa|^{1+\alpha}]=&\dfrac{1}{\sqrt{\text{det}(2\pi\Sigma)}}\times\\\
&\int_{\mathbb{R}^{5}}\left|\dfrac{x_{2}^{2}x_{3}-2x_{1}x_{2}x_{4}+x_{1}^{2}x_{5}}{(x_{1}^{2}+x_{2}^{2})^{3/2}}\right|^{1+\alpha}\exp{(-1/2(\mathbf{x}^{T}\Sigma^{-1}\mathbf{x}))}d\mathbf{x}d\mathbb{P}_{1}(\mathbf{x}^{\prime}).\end{split}$

By changing the variables to $x_{1}=r\cos{\theta},x_{2}=r\sin{\theta}$ and
keeping other variables same, we get ,

$\begin{split}\mathbb{E}[|\kappa|^{1+\alpha}]=&\dfrac{1}{\sqrt{\text{det}(2\pi\Sigma)}}\times\\\
&\int_{I}r^{-\alpha}\left|\sin^{2}{\theta}x_{3}-\sin(2\theta)x_{4}+\cos^{2}\theta
x_{5}\right|^{1+\alpha}\exp(-1/2(\mathbf{\tilde{x}}^{T}\Sigma^{-1}\mathbf{\tilde{x}}))drd\theta
d\mathbb{P}_{1}(\mathbf{x}^{\prime})\end{split}$

where $\mathbf{\tilde{x}}=(r\cos\theta,r\sin\theta)$ and
$I=[0,\infty]\times[0,2\pi]\times\mathbb{R}^{3}$. Now, for $0\leq\alpha<1$ the
above integral converges. Near the origin of $I$ convergence is taken care by
$\int_{0}^{1}r^{-\alpha}dr<\infty$ and away from origin $\exp(\cdots)$
dominates. The result follows from the fact that the vector
$(\partial_{xx}f,\partial_{xy}f,\partial_{yy}f)$ has all moments finite.

Integrability of curvature function: Consider a deterministic $C^{2}$-Morse
function $f$ on a compact domain $D\subset\mathbb{R}^{d}$. As above, at every
$x\in D$ which is a regular point of $f$, define $\kappa$ to be the divergence
of unit normal of $f$.

We prove that

$\int_{D}|\kappa|d\text{vol}<\infty.$

Note that except at critical points of $f$, $\kappa$ is continuous. So just
need to show that $\int_{B_{r}(x_{0})}|\kappa|d\text{vol}<\infty$ for a
critical point $x_{0}$ of $f$ and a small enough ball $B_{r}(x_{0})$ around
$x_{0}$.

We have $\nabla f(x)=H(f)|_{x_{0}}(x-x_{0})+O(||x-x_{0}||^{2})$, by Taylor’s
series. Since $f$ is Morse, we can invert $H(f)|_{x_{0}}$ to have

$||\nabla f(x)||\geq C\dfrac{||x-x_{0}||}{||H(f)_{x_{0}}^{-1}||}.$

Since $\partial_{xx}f,\partial_{xy}f,\partial_{yy}f$ are all bounded on $D$
and $|\partial_{x}f(x)|\leq c_{1}||x-x_{0}||$, $|\partial_{y}f(x)|\leq
c_{2}||x-x_{0}||$ near $x_{0}$ and again, exploiting the equation (14), we
have

$\int_{B_{r}(x_{0})}|\kappa|d\text{vol}\leq\tilde{C}\int_{B_{r}(x_{0})}\dfrac{1}{||x-x_{0}||}d\text{vol}.$

But we have

$\int_{B_{r}(x_{0})}\dfrac{1}{||x-x_{0}||}d\text{vol}<\infty$

for any $d\geq 2$. This completes the proof that the curvature function is
integrable on $D$.

Product of Gaussian random variables: If $X,Y$ are standard Gaussians of
correlation $\rho$, then the density of the product $Z=XY$ is given by

$\psi_{Z}(z)=\frac{1}{\pi\sqrt{1-\rho^{2}}}\exp\left(\frac{\rho
z}{1-\rho^{2}}\right)K_{0}\left(\frac{|z|}{1-\rho^{2}}\right)$ (15)

where $K_{0}$ is the Bessel function of the second kind. Integrating $\psi$
over $\mathbb{R}^{+}$,

$\mathbb{P}(Z>0)=\frac{\pi-\arccos(\rho)}{\pi}$ (16)

i.e. this probability is close to $1$ or $0$ for $\rho$ close to $\pm 1$.

## References

  * [1] J-M Azaïs and M. Wschebor, _Level sets and extrema of random processes and fields_ , John Wiley & Sons, Inc., Hoboken, NJ, 2009. MR 2478201
  * [2] J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, _The Statistics of Peaks of Gaussian Random Fields_ , The Astrophysical Journal 304 (1986), 15.
  * [3] D. Beliaev, _Gaussian fields and percolation_ , arXiv.2207.13448 (2022).
  * [4] D. Beliaev and R. W. Maffucci, _Coupling of stationary fields with application to arithmetic waves_ , Stochastic Processes and their Applications 151 (2022), 436–450.
  * [5] M. V. Berry, _Regular and irregular semiclassical wavefunctions_ , Journal of Physics A: Mathematical and General 10 (1977), no. 12, 2083–2091.
  * [6] H. B. Lawson, _Lectures on minimal submanifolds_ , 2nd ed ed., Mathematics lecture series ; 9, Publish or Perish, Berkeley, CA, 1980 (eng).
  * [7] M.S. Longuet-Higgins, _The statistical analysis of a random, moving surface_ , Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 249 (1957), 321–387.
  * [8] F. Nazarov and M. Sodin, _Asymptotic laws for the spatial distribution and the number of connected components of zero sets of gaussian random functions_ , arXiv.1507.02017 (2015).
  * [9] M. Sodin, _Lectures on random nodal portraits_ , Probability and statistical physics in St. Petersburg 91 (2016), 395–422.
  * [10] I. Wigman, _On the nodal structures of random fields – a decade of results_ , arXiv.2206.10020 (2022).
  * [11] KJ Worsley, S Marrett, P Neelin, AC Vandal, KJ Friston, and AC Evans, _A unified statistical approach for determining significant signals in images of cerebral activation_ , Human brain mapping 4 (1996), no. 1, 58—73.


# Multi-Objective Hyperparameter Tuning and Feature Selection using Filter
Ensembles

Martin Binder*
Ludwig-Maximilians-Universität München
<EMAIL_ADDRESS>&Julia Moosbauer*
Ludwig-Maximilians-Universität München
<EMAIL_ADDRESS>&Janek Thomas
Fraunhofer Institute for Integrated Circuits IIS
<EMAIL_ADDRESS>&Bernd Bischl
Ludwig-Maximilians-Universität München
<EMAIL_ADDRESS>
###### Abstract

Both feature selection and hyperparameter tuning are key tasks in machine
learning. Hyperparameter tuning is often useful to increase model performance,
while feature selection is undertaken to attain sparse models. Sparsity may
yield better model interpretability and lower cost of data acquisition, data
handling and model inference. While sparsity may have a beneficial or
detrimental effect on predictive performance, a small drop in performance may
be acceptable in return for a substantial gain in sparseness. We therefore
treat feature selection as a multi-objective optimization task. We perform
hyperparameter tuning and feature selection simultaneously because the choice
of features of a model may influence what hyperparameters perform well.

We present, benchmark, and compare two different approaches for multi-
objective joint hyperparameter optimization and feature selection: The first
uses multi-objective model-based optimization. The second is an evolutionary
NSGA-II-based wrapper approach to feature selection which incorporates
specialized sampling, mutation and recombination operators. Both methods make
use of parameterized filter ensembles.

While model-based optimization needs fewer objective evaluations to achieve
good performance, it incurs computational overhead compared to the NSGA-II, so
the preferred choice depends on the cost of evaluating a model on given data.

***These authors contributed equally to this work.

_K_ eywords Feature Selection $\cdot$ Hyperparameter Optimization $\cdot$
Multiobjective Optimization $\cdot$ Filter Ensembles $\cdot$ Evolutionary
Algorithms $\cdot$ Bayesian Optimization

##  1 Introduction

Machine learning models often need to satisfy multiple objectives
simultaneously to accomodate the nature of a practical setting. Usually the
main goal is predictive performance. Especially on large and complex datasets
this necessitates highly nonlinear algorithms, which have hyperparameters that
need to be chosen carefully. _Hyperparameter optimization_ poses a substantial
challenge in machine learning: Besides few model-specific methods [15, 29],
there are no general analytic representations of model performance w.r.t.
hyperparameter settings. Performance therefore needs to be estimated using
test-set evaluation or cross-validation. Hyperparameter optimization is
therefore an expensive black-box optimization problem.

Besides predictive performance, model sparsity is frequently another desirable
objective. According to Guyon and Elisseeff [16], sparser models help with
interpretability, i.e. a better understanding of the underlying process that
generated the data. In addition to that, predictions can be made faster and
more cost-effectively. Sparser models may even have better predictive
performance, since they regularize against overfitting.

The process of _feature selection_ aims to select a small subset of relevant
features while still constructing models with sufficient or even optimal
predictive performance. There are two distinct model-agnostic approaches to
feature selection [16]: Filters and wrappers. _Filters_ use proxy measures to
rank features by their estimated explanatory power, independently of the
learning algorithm being employed. These include information theoretic
measures, correlation measures, distance measures or consistency measures [7].
In contrast, _wrappers_ [27] optimize model test-set performance directly over
the space of feature subsets. Because every feature subset evaluation requires
either one or multiple model fits, exhaustive search is usually infeasible,
and a black-box discrete optimization search strategy is necessary. Commonly
used are simple greedy methods like forward or backward search. More advanced
methods like evolutionary algorithms can improve upon this [35]. Because they
directly optimize learner performance, wrappers often yield better results
[35].

Feature selection is often considered as a single-objective task. Sometimes
the feature selection step is only used to optimize performance [27]. However,
in many applications it is desirable to forego a small drop in performance for
a substantial gain in sparseness. This leads to a natural treatment of the
feature selection problem as a multi-objective optimization problem: Maximize
predictive performance while minimizing the number of features selected.
Feature selection methods may aggregate model performance and number of
features into a single objective function through a penalization term [35].
However, this implies a trade-off between performance and sparsity must be
specified a-priori, which may be difficult.

Multi-objective optimization methods try to find a _set_ of solutions that
represent different trade-offs between the different goals, enabling the user
to consider the possible alternatives and to choose a fitting solution
a-posteriori. This is beneficial for feature selection [34].

Hyperparameter optimization and feature selection are often performed in
separate steps. We argue that jointly optimizing over the combined spaces of
hyperparameters and feature subsets is beneficial and appropriate: The optimal
choice of hyperparameter configuration is very likely to depend on the
specific features that are included and vice-versa. Also, it is likely more
computationally efficient to explore the joint spaces simultaneously. When
using the wrapper approach this combination is not trivial: The exponentially
large binary search space of selected features now has to be fused with the
mixed numeric-categorical space of hyperparameters.

We present and adapt model-agnostic holistic approaches for both aspects
discussed above: multi-objective and joint optimization of hyperparameters and
feature sets. To guide the search, our methods make use of combinations of
feature filter scores which are combined in a _filter ensemble_.

Our approaches can be considered _hybrid_ filter-wrappers: they are
fundamentally “wrapper”-based, because they optimize model-performance, but
also make use of filters.

##  2 Related Work and Contributions

In recent literature, _Bayesian optimization_ (“BO”), also referred to as
_model-based optimization_ , has become a popular method for hyperparameter
tuning of learning algorithms [31], often outperforming simple grid search or
random search. Originally, BO generally relied on a Gaussian process (GP)
model, such as in the Efficient Global Optimization (EGO) algorithm [23].
However, the GP does not typically scale well to high dimensions and large
numbers of data points. Random forests have been proposed as an alternative
surrogate model, as used in the popular _SMAC_ hyperparameter optimization
method [21].

Automated machine learning (AutoML) deals with the configuration and
optimization of complete machine learning pipelines, often encompassing data
pre-processing, ML models, ensembling, and possibly post-processing steps.
AutoML may encompass choosing among the many different possibilities of what
method to use at each stage of the pipeline, optimizing the hyperparameters
for these methods, and combining them in ways that yield well-performing
ensembles. Feature selection by filter or wrapper methods is an important pre-
processing step already part of some AutoML frameworks. Consequently, AutoML
has the potential to jointly optimize hyperparameters and included features.
_autosklearn_ [13] for example integrates a filter-based feature selector,
parameterized by a filter measure and a percentage indicating the fraction of
highest-ranked features to be included. Bayesian optimization is used to find
the optimal filter (among a set of possible filters) and the best feature
selection rate.

However, even though hyperparameter optimization and feature selection are
both present in some AutoML frameworks, it is not their goal to find a good
trade-off between predictive performance and sparseness. Feature selection is
merely used to improve predictive performance, without considering the
preference for sparse models in light of better interpretability or other
benefits. In fact, these frameworks often tend to produce considerably complex
models. They may even introduce additional features through feature
engineering in pursuit of increasing predictive performance as the only goal.
The results are often large and heterogeneous ensembles that are hard to
interpret and deploy [30].

Evolutionary algorithms are especially well suited for multi-objective
optimization. According to the survey by Xue et al. [35], _genetic algorithms_
(GAs) are among the most commonly applied techniques for multi-objective
feature selection. Inspired by natural evolution, GAs apply recombination and
mutation operators to iteratively improve the population of solution
candidates. Through techniques like non-dominated sorting [8, 24], genetic
algorithms have become a powerful tool for multi-objective optimization.
Several GA-based methods have been proposed for the task of wrapper-based
feature selection [12, 17, 33].

There have been first investigations on simultaneous multi-objective
hyperparameter optimization and feature selection using GAs: Bouraoui et al.
[3] proposed an SVM-wrapper approach based on the NSGA-II, using a shared
representation of the feature configuration and algorithm hyperparameters.
While their method of combining the search spaces is similar to our GA-based
approach, Bouraoui et al. [3] do not use specialized initialization and
mutation operators that we have found to be necessary for good
performance111See the Ablation Study in the Supplement, where their approach,
“Variant (1)”, is outperformed by all our methods on every dataset and
learning algorithm with no exception.. Another limitation is that their
approach is not model agnostic but limited to SVMs only.

Numerous methods have been proposed for _feature filtering_ , and there are
known ways of combining filters into _filter ensembles_. These may rely on
applying the same filter on diverse datasets (“homogeneous approach”), or
different filters on the same dataset (“heterogeneous approach”), and there
are different ways of aggregating the filter rankings [2]. We extend the
heterogeneous approach into a hybrid filter-wrapper selection method by
optimizing a parameterized ranking aggregator.

There has been a lot of research on feature selection and hyperparameter
optimization. However, we found that there is no _general_ algorithm or
framework to perform model-agnostic simultaneous hyperparameter optimization
and feature selection for predictive performance and model sparsity in a
multi-objective fashion. We therefore choose to tackle this problem from two
directions. First we adapt a standard evolutionary multi-objective
optimization method, the NSGA-II, for the particular problem of feature
selection and hyperparameter tuning. We contrast this with a classical
hyperparameter optimization algorithm, based on BO, extended to perform
feature selection.

Our main contributions are:

  1. 1.

We adapt the NSGA-II to the problem at hand by introducing specialized
sampling and mutation operators that make use of filter ensembles to enhance
optimization performance.

  2. 2.

We investigate how multi-objective Bayesian optimization (MOBO) methods can be
used for this problem and propose an effective method that uses a filter
ensemble for feature selection.

  3. 3.

By conducting a benchmark of these approaches on a variety of tasks for
different machine learning algorithms, we provide a comparison between our
approaches and show how these algorithms perform compared to suitable
baselines.

##  3 Problem Statement

We consider the machine learning problem of a given feature space
$\mathcal{X}$ of vectors with $p$ components (“features”), an arbitrary
outcome space $\mathcal{Y}$ (for example $\\{-1,1\\}$ for a binary
classification task), and a performance measure
$L:\mathbb{R}^{g}\times\mathcal{Y}\rightarrow\mathbb{R}$ measuring the quality
of predictions (in $\mathbb{R}^{g}$) given ground truth values. $g$ is 1 for
regression tasks, and equal to the number possible outcome classes for
classification. Data samples
$\mathcal{D}=\left\\{\left(\mathbf{x}^{(1)},y^{(1)}\right),\ldots,\left(\mathbf{x}^{(n)},y^{(n)}\right)\right\\}\in{\left(\mathcal{X}\times\mathcal{Y}\right)}^{n}$
are assumed to be $n$ i.i.d. realizations of random variables
$(\boldsymbol{X},\boldsymbol{Y})$ which follow a joint distribution
$\mathbb{P}_{XY}$. A dataset is thus notably characterized by $n$, the number
of samples available, and $p$, the number of _features_.

Let $\mathcal{A}_{\boldsymbol{\lambda},\boldsymbol{s}}$ be a learning
algorithm that takes the given dataset $\mathcal{D}$ and constructs a model
$f:\mathcal{X}\rightarrow\mathbb{R}^{g}$. The learning algorithm is
parameterized by the _feature configuration vector_
$\boldsymbol{s}\in{\\{0,1\\}}^{p}$, where $s_{j}=1$ denotes that feature $j$
is included in the model. The sparseness of the resulting model is thus
determined by the _Hamming weight_ of $\boldsymbol{s}$, i.e. the number of
components of $\boldsymbol{s}$ that are 1. _Hyperparameters_ of the learning
algorithm are summarized in a vector $\boldsymbol{\lambda}\in\Lambda$.
$\Lambda$ is a (possibly mixed) bounded space and may contain numeric,
integer, and categorical values222Hierarchical dependencies of
hyperparameters, e.g. kernel hyperparameters that appear only if a specific
kernel is chosen, are a common extension but not considered in this work..

In general we are trying to construct a model
$f=\mathcal{A}_{\boldsymbol{\lambda},\boldsymbol{s}}(\mathcal{D})$ which
minimizes the _generalization error_
$\textrm{GE}[f]=\mathbb{E}_{\mathbb{P}_{XY}}\left[L(f(\mathbf{x}),y)\right]$.
However, the generalization error can only be estimated using in-sample data
$\widehat{\textrm{GE}}[\mathcal{A}_{\boldsymbol{\lambda},\boldsymbol{s}},\mathcal{D}]$
through a resampling technique such as cross-validation.

This setup gives rise to the bi-objective hyperparameter optimization and
feature selection problem:

$\min_{\boldsymbol{\lambda}\in\Lambda,\boldsymbol{s}\in{\\{0,1\\}}^{p}}\left(\widehat{\textrm{GE}}\left[\mathcal{A}_{\boldsymbol{\lambda},\boldsymbol{s}},\mathcal{D}\right],\sum_{i=1}^{p}c_{i}s_{i},\right)\textrm{.}$

The setup regards estimated generalization error as one objective and the cost
of features considered as another.

It is possible, and a trivial extension of our method, to consider arbitrary
costs $c_{i}$ for each feature $i$. However, in our benchmarks we limit
ourselves to equal costs $c_{i}=1/p$. The resulting measure corresponds to the
_fraction of selected features_ , $\mathrm{ffrac}=\frac{1}{p}\sum_{i}s_{i}$,
ranging from 0 to 1.

##  4 Multi-Objective Hyperparameter Tuning and Feature Selection: Two
Approaches

Our two methods are based on two popular multi-objective optimization methods,
adapted for the particular task of hyperparameter tuning and feature
selection: A model-based approach, and an approach based on an evolutionary
algorithm. Both approaches make use of feature filters to accelerate the
search.

###  4.1 Feature Filters and Filter Ensembles

The number of possible feature configurations $\boldsymbol{s}$ is exponential
in $p$, so a large number of performance evaluations would be necessary to
explore the performance space. Therefore, both our optimization approaches
make use of _feature filters_. These are methods that heuristically score
individual features according to their apparent relevance for the outcome
variable. Because there are various methods to estimate this relevance, we
consider a collection of $M$ feature filters, which generate scores
${F^{m}(\mathcal{D})}_{j}$, $m=1,\ldots,M$, $j=1,\ldots,p$ for each feature
$j$ of a training dataset $\mathcal{D}$. Different methods may generate scores
on different scales, but we rank-transform and scale these scores to values
ranging in equi-distant steps from least (value of 0) to most (value of 1)
relevant for the outcome variable.

We propose combining multiple filter methods into _filter ensembles_ similarly
to Dittman et al. [11], but extending their method by using _weighted_ average
rank aggregation. The ensemble filter score $\mathrm{EF}_{j}$ for feature $j$
is calculated as the weighted average according to weight vector
$\boldsymbol{w}\in[0,1]^{M},\sum_{i}w_{i}=1$:

$\mathrm{EF}_{j}(\boldsymbol{w})=\sum_{m=1}^{M}w_{m}{F^{m}(\mathcal{D})}_{j}\textrm{.}$
(1)

The weighting parameter $\boldsymbol{w}$ can be optimized, extending the
filter into a hybrid filter-wrapper approach.

###  4.2 Bayesian Optimization Approach

_Bayesian optimization_ (BO), also often referred to as _sequential model-
based optimization_ , has been successfully used for machine learning
hyperparameter optimization in many applications [31]. The principle of BO is
based on two steps, which are performed in turn. First, a so-called _surrogate
model_ is fitted to model the relationship between decision variables (e.g.
hyperparameter values) and the objective value (e.g. estimated generalization
performance). The surrogate model generates cheap approximations of the
(generally expensive to evaluate) objective function values. In a second step,
an _infill criterion_ is used to find promising decision values to be
evaluated on this expensive function. The infill criterion has to face a
trade-off between “exploitation”—evaluating points for which the surrogate
predicts good performance—and “exploration”—evaluating points where predictive
uncertainty is high.

There are different ways of adapting BO to perform multi-objective Bayesian
optimization (“MOBO”) [20], which has been applied successfully to
hyperparameter tuning [18]. We chose the _Parego_ method [26]. Parego
scalarizes the different objectives through the Chebyshev norm using a random
weight vector, which is sampled again for every point being proposed. Parego
has the advantage over many other MOBO methods of only requiring a single
surrogate model fit for a proposed point [18], but other MOBO approaches can
trivially be used as a substitute within our method.

The number of possible feature configurations is exponential in $p$, so a
large number of evaluations would be necessary to fit an accurate surrogate
model on performance values. We therefore investigate the use of prior
knowledge from _feature filter_ scores to reduce the dimensionality of the
search space. Two possible methods for simultaneous use of multiple feature
filter methods are considered (see also Figure 1):

GA-MO-FE: NSGA-II with feature ensemble mutation

GA-MO: NSGA-II with Hamming-weight preserving mutation

GA-MO-FE-NJ: NSGA-II for features; hyperparameters fixed

BO-MO-FE, BO-SO-FE, BO-MO-FE-NJ: BO with filter ensemble

BO-MO, BO-SO: BO with individual filter selection

Figure 1: Representation of individuals in different variants of our proposed
optimization methods as described in Section 5.3; $\mathbf{\lambda}$,
$\mathbf{w}$, and $\mathbf{s}$ are vectors as described in Section 4.
Figure 2: Example population (800 samples) of _xgboost_ configurations
randomly sampled and evaluated on the _madelon_ task. Red:
$\frac{1}{2}$-Bernoulli sampling. Blue: Geometric initialization (Section
4.3.1). Triangles show non-dominated individuals. Although randomly
initialized and not yet optimized, the blue individuals dominate a larger
area.

####  4.2.1 Individual filter selection

In this method, we introduce a discrete _filter selection hyperparameter_ $m$,
as well as a _feature fraction hyperparameter_ $\mathrm{ffrac}\in[0,1]$. For
each model evaluation, only the most relevant $\lceil
p\cdot{}\mathrm{ffrac}\rceil$ features, according to filter with index $m$,
are included in the model. This approach is similar to the one taken in auto-
sklearn [13], although there the feature selection problem was not considered
as a multi-objective problem.

####  4.2.2 Filter ensemble selection

This method uses the filter ensemble as shown in Equation 1 and introduces the
vector $\boldsymbol{w}$, as well as the aforementioned feature fraction
$\mathrm{ffrac}$, as hyperparameters. The most relevant $\lceil
p\cdot{}\mathrm{ffrac}\rceil$ features, according to
$\mathrm{EF}_{j}(\boldsymbol{w})$, are included in the model.

###  4.3 Evolutionary Approach

The _Nondominated Sorting Genetic Algorithm II_ (NSGA-II) [9] is an
evolutionary multi-objective algorithm that uses nondominated-sorting to
preferably select individuals close to the Pareto-front of the problem. It
iterates through generations of each a _reproduction_ , a _crossover_ , a
_mutation_ , and a _survival_ step that generate the population of the next
generation.

GAs often represent individuals as vectors of binary, discrete, or continuous
values, depending on the optimization problem. Because hyperparameters
generally have various types, we use the Cartesian products of operators that
operate in different ways on the various types, following Li et al. [28]. This
means e.g. that numeric hyperparameters undergo Gaussian mutation, while
categorical hyperparameters undergo uniform mutation etc. Table 1 summarizes
the chosen recombination and mutation operators for the respective
hyperparameter types. For the hyperparameter mutations, we use self-adapting
step sizes and mutation probabilities as suggested by Li et al. [28].

The feature configuration vector parameter $\mathbf{s}$ plays a special role,
because $\mathbf{s}$ maps to the objective of the fraction of selected
features, which is being optimized, in a straightforward manner. For datasets
with many features it makes up a large part of the search space under
consideration compared to the other hyperparameters. The initialization and
mutation performed on this parameter should therefore be considered in detail.

Table 1: Summary of mutation and recombination parameters used in our GA-based algorithms. _s.a._ : the parameter is controlled by self-adaption [28]. Filter-ensemble mutation of feature configurations is used in GA-MO-FE(-NJ), Hamming-weight preserving mutation in GA-MO. parameter | recombination |  mutation
---|---|---
numeric | SBX $(\eta=5)$ |  Gaussian
($p=0.1$, $\sigma^{2}$ = s.a.)
integer | rounded SBX $(\eta=5)$ |  rounded Gaussian
($p=0.1$, $\sigma^{2}$ = s.a.)
categorical | uniform ($p=0.5$) |  uniform mut. ($p$ = s.a.)
binary | uniform ($p=0.5$) |  uniform mut. ($p$ = s.a.)
features | uniform ($p=0.5$) |  Filter-ensemble or
Hamming-weight pres.

####  4.3.1 Geometric initialization

A naive approach for feature configuration initialization would be Bernoulli-
sampling of each feature selection bit $s_{i}$ individually, possibly biased
towards a low expected number of selected features to favor relatively sparse
solutions [1]. However, this gives rise to a binomial distribution of the
number of selected features $\sum_{i}s_{i}$ with standard deviation $\sim
O\left(\sqrt{p}\right)$. For even moderately large values of $p$, this fails
to cover the objective space evenly along the dimension of the selected
feature fraction, see Figure 2. We propose to sample values of
$\boldsymbol{s}$ such that the sum of selected features covers the whole
feature fraction objective. Therefore, we elect to use a truncated geometric
distribution of number of included features to encode our preference for
sparse models. This is achieved by sampling the desired number of included
features $S$ as a truncated geometrically distributed random integer between
$0$ and $p$, and then uniformly sampling from all vectors $\mathbf{s}$ that
satisfy $\sum_{i}s_{i}=S$.

The method introduces the success probability of the geometric distribution as
a configuration parameter. It can be set by the user to encode a relative
preference for sparsity. We chose to use an empirically determined value by
fitting decision trees on 100 random subsets of 90% of the data set and
determining the average number of distinct split variables.

####  4.3.2 Filter-ensemble based initialization

It is possible to enhance the geometric initialization by including prior
knowledge gained from feature filter methods. The goal is to select the
features that have high filter scores with larger probability than the ones
with lower ranking, while still having an approximately geometric distribution
over the number of total features selected. We make the initial distribution
of bit $j$ dependent on the filter ensemble value as described in equation 1,
with $\boldsymbol{w}$ uniformly randomly sampled from the simplex
$\boldsymbol{w}\in{\left[0,1\right]}^{M}$, $\sum_{m=1}^{M}w_{m}=1$. For each
individual to be initialized, we first sample $S$ as in 4.3.1. Each bit $j$ is
then sampled from a Bernoulli-distribution with parameter:

$\pi_{\mathrm{B}}(\boldsymbol{w},S)=\frac{\mathrm{EF}_{j}(\boldsymbol{w})\left(S+1\right)}{\mathrm{EF}_{j}(\boldsymbol{w})S+\left(1-\mathrm{EF}_{j}(\boldsymbol{w})\right)\left(p-S\right)+1}\textrm{.}$
(2)

####  4.3.3 Hamming-weight preserving mutation

Performing random bit-flip mutation on the feature selection vector entails
similar problems to Bernoulli-initialization: a bit-flip with probability
$\pi_{\mathrm{Mut}_{\mathrm{bin}}}$ is equivalent to erasing a bit with
probability $2\pi_{\mathrm{Mut}_{\mathrm{bin}}}$ and sampling it anew from a
$\frac{1}{2}$-Bernoulli distribution. This biases the mutation result towards
the $\frac{1}{2}$ point of the feature fraction objective. Instead, we choose
to approximately preserve the Hamming weight. This is done by sampling erased
bits from a Bernoulli-distribution with parameter
$\pi_{\mathrm{B}}=(S+1)/(p+2)$, where $S=\sum_{i}s_{i}$ is the Hamming weight
of the original vector.

####  4.3.4 Filter-ensemble based mutation

Just as for initialization, mutation can be made dependent on feature filter
ensemble ranks to preferentially include more relevant features. As in the
Hamming-weight preserving mutation, each bit is erased with a probability of
$2\pi_{\mathrm{Mut}_{\mathrm{bin}}}$ and then drawn from a Bernoulli-
distribution. The parameter is as in equation 2 for filter-based
initialization, with $S$ the Hamming weight of the original vector, and
$\boldsymbol{w}$ a weighting vector. This weighting vector itself becomes part
of the search space to achieve self-adaption, as described in Li et al. [28],
which results in $\mathbf{w}$ being optimized during the NSGA-II run. This
feature configuration mutation method is (approximately) Hamming-weight
preserving.

###  4.4 Implementation and Reproducibility

All proposed methods are implemented and publicly accessible through the
mosmafs R package published on
CRAN333https://CRAN.R-project.org/package=mosmafs. For full reproducibility of
the results we publish the code used to perform the benchmark experiments on
Github444https://github.com/compstat-lmu/mosmafs/tree/master.

##  5 Experiments

We conduct experiments to answer the following research questions empirically:

  1. 1.

_Evolutionary vs. Bayesian Optimization_ : How do the proposed methods—the
evolutionary and the Bayesian optimization approach—perform relative to each
other?

  2. 2.

_Effect of Filter Ensembles_ : Do the methods benefit from using _filter
ensembles_?

  3. 3.

_Multi-Objective vs. Single-Objective_ : Does multi-objective optimization
find much sparser solutions without a major loss in predictive performance
compared to single-objective optimization?

  4. 4.

_Simultaneous Hyperparameter Tuning and Feature Selection_ : Is it beneficial
to perform hyperparameter optimization and feature selection _simultaneously_
compared to performing the tasks sequentially?

###  5.1 Benchmark Datasets

We consider real-world binary classification tasks that are publicly
accessible through the OpenML platform [32] (see Table 2). To eliminate
algorithmic factors that might influence the result (like class imbalance
correction, feature encoding, or handling of missing values), we included
datasets fulfilling the following criteria: (roughly) balanced outcome
classes, a purely numeric feature space, and no missing values. Further,
datasets have been chosen to represent a large variety in terms of
dimensionality ($30\leq p\leq 10937$)555To limit the computational resources
necessary for our experiments we did not investigate even larger datasets. and
in terms of instances per dimension ($0.02\leq n/p\leq 18.97$).

Although our method can be used for regression as well as classification
tasks, we did not perform experiments on the former.

Table 2: Description of the datasets being used. $n$ denotes the number of observations and $p$ the total number of features. The class ratio gives the proportional size of the smaller outcome class. The dataset id (did) is the unique identifier for the dataset on the OpenML platform [32].  name | $n$ | $p$ | class ratio | $n/p$ | did
---|---|---|---|---|---
wdbc | $569$ | $30$ | $0.37$ | $18.97$ | $1510$
ionosphere | $351$ | $33$ | $0.36$ | $10.64$ | $59$
sonar | $208$ | $60$ | $0.47$ | $3.47$ | $40$
hill-valley | $1212$ | $100$ | $0.50$ | $12.12$ | $1479$
tecator | $240$ | $124$ | $0.43$ | $1.94$ | $851$
semeion | $319$ | $256$ | $0.50$ | $1.25$ | $41973$
madeline | $3140$ | $259$ | $0.50$ | $12.12$ | $41144$
lsvt | $126$ | $307$ | $0.33$ | $0.41$ | $1484$
madelon | $2600$ | $500$ | $0.50$ | $5.20$ | $1485$
isolet | $600$ | $617$ | $0.50$ | $0.97$ | $41966$
cnae-9 | $240$ | $282$ | $0.50$ | $0.85$ | $41967$
arcene | $200$ | $9961$ | $0.44$ | $0.02$ | $1458$
AP_Breast_Colon | $630$ | $10935$ | $0.45$ | $0.06$ | $1145$
AP_Colon_Kidney | $546$ | $10935$ | $0.48$ | $0.05$ | $1137$

###  5.2 Learning Algorithms and their Hyperparameters

We consider three different classifiers that are tuned by the optimization
algorithms proposed in Section 4: The support vector machine classifier (SVM)
[5] with Gaussian kernel, extreme gradient boosting (xgboost) [4] and the
kernelized k-nearest-neighbor classifier (kknn) [37]. We decided to consider
these learners because they represent three very distinct learning paradigms,
and because they are generally regarded responsive to hyperparameter tuning.
The hyperparameter spaces that are being tuned over are presented in Table 3.
They also react differently to high dimensionality: kknn is vulnerable to the
“curse of dimensionality”, the SVM is a regularized modeling algorithm, and
the xgboost algorithm does implicit feature selection because it is a tree-
based learner.

Table 3: Hyperparameter spaces over which tuning was performed for the three learning algorithms _xgboost_ , _support vector machine_ and _kernelized k-Nearest-Neighors_ respectively. Ranges marked with ∗ were tuned on a logarithmic scale. For the large datasets madeline, madelon, arcene, AP_Breast_Colon and AP_Colon_Kidney, the xgboost _nrounds_ parameter was set fixed to 2000 and early stopping after 10 rounds was enabled.  SVM |  | xgboost
---|---|---
kernel | rbfdot |  | nrounds | $\\{1,2,\ldots,2000\\}$
sigma | $[2^{-10},2^{10}]~{}^{*}$ |  | eta | $[0.01,0.2]$
C | $[2^{-10},2^{10}]~{}^{*}$ |  | gamma | $[2^{-7},2^{6}]~{}^{*}$
|  |  | max_depth | $\\{3,\ldots,20\\}$
kknn |  | colsample_bytree | $[0.5,1]$
k | $\\{1,2,\ldots,50\\}$ |  | colsample_bylevel | $[0.5,1]$
distance | $[1,100]$ |  | lambda | $[2^{-10},2^{10}]~{}^{*}$
kernel | {rectangular, |  | alpha | $[2^{-10},2^{10}]~{}^{*}$
| optimal, |  | subsample | $[0.5,1]$
| triangular, |  |  |
| biweight} |  |  |

###  5.3 Algorithms

We perform hyperparameter optimization and feature selection for the learners
described in Section 5.2. The datasets used are presented in Section 5.1. We
use different configurations of our methods for comparisons described in the
following. See also Figure 1 for a schematic representation of individual
configuration vectors used in each method.

To answer research question 1, comparing the evolutionary and Bayesian
optimization methods, we consider the following algorithms:

GA-MO-FE

The NSGA2 with filter ensemble based initialization (as in 4.3.2) _and_
mutation (as in 4.3.4).

BO-MO-FE

Multi-objective Bayesian optimization with filter ensemble selection as
described in 4.2.2.

To study research question 2, i.e. the contribution of filter ensembles to the
methods, we compare both methods to these algorithms without filter ensembles:

GA-MO

The NSGA2 with filter ensemble based initialization (as in 4.3.2) but
_without_ filter-based mutation, instead only using Hamming-weight preserving
mutation as in 4.3.3.

BO-MO

Multi-objective Bayesian optimization with individual filter selection, as
described in 4.2.1.

We illuminate research question 3 by considering feature selection and
hyperparameter optimization as a single-objective task, run with the following
methods:

BO-SO

Single-objective Bayesian optimization with individual filter selection (see
4.2.1)—this is a method similar to auto-sklearn’s approach [13] which we
consider to be a state-of-the-art approach.

BO-SO-FE

Single-objective Bayesian optimization with filter ensemble selection (see
4.2.2), which is our strongest single objective baseline.

To address research question 4, we consider variants of our algorithms that
perform feature selection and hyperparameter optimization in separate steps.
We construct two algorithms to approximate the most straightforward
simplifications to non-joint optimization. Note that the GA-approach has a
focus on feature selection, while the BO-approach is more directed at
hyperparameter tuning.

GA-MO-FE-NJ

We initially optimize hyperparameters for the learning algorithms on each
dataset by running standard single objective Bayesian optimization on the
_full_ datasets without performing feature selection, through 500 model
evaluations. The resulting hyperparameter values are then fixed, while our
full NSGA2-variant (using both filter ensemble initialization and mutation)
only performs multi-objective feature selection.

BO-MO-FE-NJ

To optimize hyperparameters and filter weights we use the BO-SO-FE method, and
then evaluate the final model with the fixed hyperparameters. We iterate over
equally spaced $\mathrm{ffrac}$ values and evaluate the model with all
filters. This constructs a set of models with different trade-offs between
sparsity and model performance, from which we construct a Pareto-set.

The methods have access to $M=5$ filter methods666The methods were chosen by
running many filters on a range of datasets and performing hierarchical
clustering on the differences of their feature rankings; see the Supplement
for more details.: Information Gain, Random Forest Feature Importance [22],
Joint Mutual Information (JMI) [36], Minimal Conditional Mutual Information
Maximization (CMIM) [14], Area Under the Curve (AUC).

All NSGA2 variants use $\mu=80$ as population and $\lambda=15$ as offspring
size which Khan and Baig [25] found to perform well in a feature selection
setting. We choose an overall per-individual mutation probability of 0.3 and
an overall per-pair-of-individuals crossover probability of 0.7.

Our Bayesian optimization methods use a random forest as surrogate to model
the mixed discrete and continuous hyperparameter space similar to the state-
of-the-art hyperparameter optimization toolbox SMAC [21]. The infill criterion
used is LCB [6]. In each iteration, a batch of $15$ configurations is proposed
in parallel as described in Horn et al. [19].

###  5.4 Evaluation

We measure the performance of resulting models by their mean misclassification
error
($\mathrm{mmce}=\frac{1}{n}\sum_{i}\mathbb{I}_{y^{(i)}\neq\hat{y}^{(i)}}$) on
a validation set with ground truth values $y^{(i)}$ and model predictions
$\hat{y}^{(i)}$.

To get an estimate of the optimization performance that is unbiased by
potential overtuning, we performed _nested resampling_ : During the whole
optimization run, each optimization algorithm is only allowed to assess model
performance $\mathrm{mmce}_{\mathrm{optim}}$ through (_inner_) cross-
validation on a optimization set $\mathcal{D}_{\mathrm{optim}}$. The final
performance of the algorithm is reported as the performance of the solution
candidates trained on $\mathcal{D}_{\mathrm{optim}}$ and evaluated on a test
set $\mathcal{D}_{\mathrm{test}}$. This procedure is repeated $10$ times on
_outer_ cross-validation folds
$\left(\mathcal{D}^{(k)}_{\mathrm{optim}},\mathcal{D}^{(k)}_{\mathrm{test}}\right),k=1,...,10$.

As a proper multi-objective performance measure we consider the _dominated
hypervolume_ [38] (domHV) with reference point $W=(1,1)$, which corresponds to
the worst possible values w.r.t. the two objectives. To prevent overtuning
effects from skewing our results, we proceed as follows:

  1. 1.

We let each optimization algorithm report its Pareto set, taking into account
only $\mathcal{D}_{\mathrm{optim}}$. This corresponds to returning non-
dominated individuals w.r.t.
$(\mathrm{mmce}_{\mathrm{optim}},\mathrm{ffrac})$. GA-based methods take only
candidates from their current generation, while BO-based methods get the
Pareto set of all candidates seen so far.

  2. 2.

We calculate the _generalization_ dominated hypervolume
$\textrm{domHV}_{\mathrm{gen}}$, i.e. the hypervolume that is dominated by
this Pareto set of candidate configurations w.r.t. the
$(\mathrm{mmce}_{\mathrm{test}},\mathrm{ffrac})$ evaluated on
$\mathcal{D}_{\mathrm{test}}$.777Note that $\textrm{domHV}_{\mathrm{gen}}$ is
not the same as the dominated hypervolume of a population on the test set.
Instead, only the individuals that are non-dominated according to
$\mathcal{D}_{\mathrm{optim}}$ are used to calculate their dominated
hypervolume on $\mathcal{D}_{\mathrm{test}}$.

For each of the experiments, we allow a maximum number of $2000$ model
evaluations, with one evaluation corresponding to computing a single inner
$10$-fold cross-validation on $\mathcal{D}_{\mathrm{optim}}$. This corresponds
to approximately $130$ iterations for both the model-based and the
evolutionary approach.

##  6 Results

###  6.1 Evolutionary vs. Bayesian Optimization, Effect of Filter Ensembles

The global rank analysis as well as a critical difference test [10] presented
in Figure 3 show that our most advanced BO-based method, BO-MO-FE,
significantly outperforms all GA-based methods. The optimization trace seems
to indicate that including the feature ensemble also confers an advantage,
although this is not statistically significant after 2000 evaluations.

Figure 4 shows an additional aspect to take into consideration: overall
runtime. BO-MO-FE usually has a slight performance advantage over GA-MO-FE in
absolute terms, but comes with computational overhead that may have to be
considered if individual model performance evaluations are relatively fast.
The overhead becomes less relevant when optimizing large datasets or slow
models.

###  6.2 Multi-Objective vs. Single-Objective

The goal of our multi-objective methods is not to optimize predictive
performance by itself, and instead to explore the possible trade-offs between
performance and sparseness. It is still interesting to look at the best
performing model configurations being found by the multi-objective methods,
and to compare them to the best models found by single-objective methods, as
done in Figure 5. Here, again, the BO-based methods outperform GA-based
methods. Except for the _kknn_ learning algorithm, the MO-method is on par
with the SO methods.

The BO-SO method employs a feature filtering step, and the models it chooses
will often be sparse to some degree, for regularization. Figure 6 compares the
BO-SO method to the multi-objective methods in terms of both sparseness and
predictive performance: For each learning algorithm, dataset, and cross-
validation fold of BO-SO, the best performing model from BO-MO-FE and GA-MO-FE
were chosen that are _at least_ as sparse. The plot shows that both MO methods
often found much sparser models than the SO baseline while giving up very
little in predictive performance.

###  6.3 Simultaneous Hyperparameter Tuning and Feature Selection

Table 4 shows the $\textrm{domHV}_{\mathrm{gen}}$ of our joint optimization
methods, compared to their non-joint correspondents. In most cases the joint
optimization confers a considerable performance advantage, especially compared
to the BO-SO-FE-NJ method. We assume that the advantage of joint over non-
joint methods depends on how much the hyperparameter performance of a model
interacts with sparseness.

Figure 3: Comparison of different multi-objective optimization methods. Top:
the results of a global rank analysis based on $\textrm{domHV}_{\mathrm{gen}}$
to compare the proposed methods BO-MO(-FE) and GA-MO(-FE), i.e. both the BO
and GA with and without filter ensemble. Ranks are computed per dataset and
algorithm (ties are ranked by their average rank) and then averaged. Higher
values are better. Bottom: Non-parametric critical difference test [10]
performed after the full budget of 2000 evaluations at significance level
$\alpha=0.05$, based on $\textrm{domHV}_{\mathrm{gen}}$ rank as above. The BO-
MO-FE method statistically significantly outperforms all other shown methods,
although the absolute difference is small (see Figure 4).
Figure 4: Performance ($\textrm{domHV}_{\mathrm{gen}}$) and runtime of GA-MO-
FE (tail end of arrow) and BO-MO-FE (head of arrow) on each dataset (Table 2)
and learning algorithm (Table 3), averaged over 10 outer cross-validation
runs. BO-MO-FE has moderately improved performance over GA-MO-FE, but at a
cost in runtime.
Figure 5: Comparison of multi-objective and single-objective methods. Shown is
the $\mathrm{mmce}_{\mathrm{test}}$ of the individual with best
$\mathrm{mmce}_{\mathrm{optim}}$ found until a particular evaluation, averaged
over 10 outer cross-validation folds and over all datasets. Lower values are
better. Note that $\mathrm{mmce}_{\mathrm{test}}$ is not monotonically
improving because individuals with better $\mathrm{mmce}_{\mathrm{optim}}$ may
be found that perform worse on the test set than previous individuals.
Figure 6: Comparison of multi-objective and single-objective methods applied to the SVM learning algorithm: Performance $\mathrm{mmce}_{\mathrm{test}}$ and the fraction of included features $\mathrm{ffrac}$ found after 2000 evaluations by baseline BO-SO (tail end of arrows), BO-MO-FE (head of solid arrows) and GA-MO-FE (head of dashed arrows). Each dataset (Table 2) is shown, values are averaged over 10 outer CV runs. Choice of individuals for MO methods described in Section 6.2. Table 4: $\textrm{domHV}_{\mathrm{gen}}$ after $2000$ evaluations of simultaneous hyperparameter tuning / feature selection methods compared to corresponding non-simultaneous methods. Best results for GA and BO in bold. Values averaged over datasets, see supplement for results by dataset.  Learner | BO-MO-FE | BO-MO-FE-NJ | GA-MO-FE | GA-MO-FE-NJ
---|---|---|---|---
kknn | 0.9217 | 0.8670 | 0.9108 | 0.9134
SVM | 0.9331 | 0.8719 | 0.9241 | 0.8892
xgboost | 0.9164 | 0.8940 | 0.9165 | 0.9121

##  7 Conclusion

Our results show that both evolutionary approaches and model-based approaches
can efficiently perform model-agnostic multi-objective optimization to
simultaneously tune hyperparameters and select features. We have also shown
that performing these tasks simultaneously has an advantage over running them
separately.

Using model-based optimization to tune over the type of filter measures as
well as fraction of included features is used in some AutoML-frameworks [13].
Our experiments suggest that the adaption of this to multi-objective
optimization works well, although using a parameterized filter ensemble
outperforms it.

Our other proposed method based on an NSGA-II enhanced with specialized
initialization and mutation sampling seems to perform almost as well, although
it does not reach the level of the best Bayesian optimization based approach.
The advantage of the NSGA-II, however, is that it does not introduce as much
computational overhead, which can make up a considerable part of overall
runtime if model performance evaluations themselves are cheap.

Our final recommendation is therefore to use the Bayesian optimization
approach in combination with parameterized Filter ensembles if model
evaluations are expensive, if computational resources are cheap, and if it is
important to get configurations that perform very close to optimal, given
their sparseness. The NSGA-II approach is suitable if model evaluations are
cheap and if marginal degradation of performance are acceptable.

## Acknowledgments

This work has been partially funded by the German Federal Ministry of
Education and Research (BMBF) under Grant No. 01IS18036A. The authors of this
work take full responsibilities for its content.

This work was supported by the Bavarian Ministry for Economic Affairs,
Infrastructure, Transport and Technology through the Center for Analytics-
Data-Applications (ADA-Center) within the framework of “BAYERN DIGITAL II”.

## References

  * Bischl et al. [2010] B. Bischl, I. Vatolkin, and M. Preuss.  Selecting small audio feature sets in music classification by means of asymmetric mutation.  In _Proceedings of the 11th International Conference on Parallel Problem Solving from Nature_ , volume 6238, pages 314–323, 2010.  ISBN 3642158439.
  * Bolón-Canedo [2018] A. Bolón-Canedo, Verónica; Alonso-Betanzos.  Ensembles for feature selection: A review and future trends.  _Information Fusion_ , 11 2018.  doi: 10.1016/j.inffus.2018.11.008.
  * Bouraoui et al. [2018] A. Bouraoui, S. Jamoussi, and Y. BenAyed.  A multi-objective genetic algorithm for simultaneous model and feature selection for support vector machines.  _Artificial Intelligence Review_ , 50(2):261–281, 2018.  ISSN 15737462.
  * Chen and Guestrin [2016] T. Chen and C. Guestrin.  _XGBoost: A Scalable Tree Boosting System_.  ACM, 2016.  ISBN 9781450342322.
  * Cortes and Vapnik [1995] C. Cortes and V. Vapnik.  Support-Vector Networks.  _Machine Learning_ , 20(3):273–297, 1995.  ISSN 15730565.
  * Cox and John [1992] D. D. Cox and S. John.  A statistical method for global optimization.  In _[Proceedings] 1992 IEEE International Conference on Systems, Man, and Cybernetics_ , pages 1241–1246. IEEE, 1992.
  * Dash and Liu [1997] M. Dash and H. Liu.  Feature selection for classification.  _Intelligent Data Analysis_ , 1(3):131–156, 1997\.  ISSN 1088467X.
  * Deb and Agrawal [1994] K. Deb and R. B. Agrawal.  Simulated binary crossover for continuous search space.  _Complex Systems_ , 9(2):115–148, 1994.  ISSN 08912513.
  * Deb et al. [2002] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan.  A fast and elitist multiobjective genetic algorithm: NSGA-II.  _IEEE Transactions on Evolutionary Computation_ , 6(2):182–197, 2002.  ISSN 1089778X.
  * Demšar [2006] J. Demšar.  Statistical comparisons of classifiers over multiple data sets.  _J. Mach. Learn. Res._ , 7:1–30, Dec. 2006.  ISSN 1532-4435.
  * Dittman et al. [2012] D. J. Dittman, T. M. Khoshgoftaar, R. Wald, and A. Napolitano.  Comparing two new gene selection ensemble approaches with the commonly-used approach.  In _2012 11th International Conference on Machine Learning and Applications_. IEEE, Dec. 2012.  doi: 10.1109/icmla.2012.175.  URL https://doi.org/10.1109/icmla.2012.175.
  * Emmanouilidis et al. [2000] C. Emmanouilidis, A. Hunter, and J. Macintyre.  A multiobjective evolutionary setting for feature selection and a commonality-based crossover operator.  In _Proceedings of the 2000 Congress on Evolutionary Computation. CEC00 (Cat. No.00TH8512)_ , volume 1, pages 309–316. IEEE, 2000.  ISBN 0-7803-6375-2.
  * Feurer et al. [2015] M. Feurer, A. Klein, K. Eggensperger, J. Springenberg, M. Blum, and F. Hutter.  Efficient and Robust Automated Machine Learning, 2015.  URL http://papers.nips.cc/paper/5872-efficient-and-robust-automated-machine-learning.
  * Fleuret [2004] F. Fleuret.  Fast Binary Feature Selection with Conditional Mutual Information.  _Journal of Machine Learning Research_ , 5:1531–1555, 2004\.
  * Franceschi et al. [2018] L. Franceschi, P. Frasconi, S. Salzo, R. Grazzi, and M. Pontil.  Bilevel programming for hyperparameter optimization and meta-learning.  In _Proceedings of the 35th International Conference on Machine Learning (ICML 2018_ , Proceedings of Machine Learning Research. PMLR, 2018.
  * Guyon and Elisseeff [2003] I. Guyon and A. Elisseeff.  An Introduction to Variable and Feature Selection.  _Journal of Machine Learning Research_ , 3:1157–1182, 2003\.  ISSN 00032670.
  * Hamdani et al. [2007] T. M. Hamdani, J.-M. Won, A. M. Alimi, and F. Karray.  Multi-objective Feature Selection with NSGA II.  In _Adaptive and Natural Computing Algorithms_ , pages 240–247. Springer Berlin Heidelberg, Berlin, Heidelberg, 2007.
  * Horn and Bischl [2016] D. Horn and B. Bischl.  Multi-objective parameter configuration of machine learning algorithms using model-based optimization.  In _2016 IEEE Symposium Series on Computational Intelligence (SSCI)_ , pages 1–8. IEEE, 2016.
  * Horn et al. [2015a] D. Horn, T. Wagner, D. Biermann, C. Weihs, and B. Bischl.  Model-based multi-objective optimization: Taxonomy, multi-point proposal, toolbox and benchmark.  In _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ , volume 9018, pages 64–78, 2015a.  ISBN 9783319159331.
  * Horn et al. [2015b] D. Horn, T. Wagner, D. Biermann, C. Weihs, and B. Bischl.  Model-based multi-objective optimization: Taxonomy, multi-point proposal, toolbox and benchmark.  In _Lecture Notes in Computer Science_ , pages 64–78. Springer International Publishing, 2015b.  doi: 10.1007/978-3-319-15934-8_5.  URL https://doi.org/10.1007/978-3-319-15934-8_5.
  * Hutter et al. [2011] F. Hutter, H. H. Hoos, and K. Leyton-Brown.  Sequential Model-Based Optimization for General Algorithm Configuration.  In _International conference on learning and intelligent optimization_. Springer, 2011.  URL http://www.cs.ubc.ca/labs/beta/Projects/SMAC/papers/11-LION5-SMAC.pdf.
  * Ishwaran et al. [2010] H. Ishwaran, U. B. Kogalur, E. Z. Gorodeski, A. J. Minn, and M. S. Lauer.  High-dimensional variable selection for survival data.  _Journal of the American Statistical Association_ , 105(489):205–217, 2010.  ISSN 01621459.
  * Jones et al. [1998] D. R. Jones, M. Schonlau, and W. J. Welch.  Efficient Global Optimization of Expensive Black-Box Functions.  _Journal of Global Optimization_ , 13:455–492, 1998.  URL http://www.ressources-actuarielles.net/EXT/ISFA/1226.nsf/0/f84f7ac703bf5862c12576d8002f5259/{$}FILE/Jones98.pdf.
  * K. Deb, A. Pratap, S. Agarwal [2002] T. M. K. Deb, A. Pratap, S. Agarwal.  A fast and elitist multi-objective genetic algorithm: NSGA-II, IEEE Transactions on Evolutionary Computation.  _IEEE Transactions on Evolutionary Computation_ , 6(2):182–197, 2002.
  * Khan and Baig [2015] A. Khan and A. R. Baig.  Multi-objective feature subset selection using non-dominated sorting genetic algorithm.  _Journal of Applied Research and Technology_ , 13(1):145–159, 2015.  ISSN 16656423.
  * Knowles [2006] J. Knowles.  ParEGO: a hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems.  _IEEE Transactions on Evolutionary Computation_ , 10(1):50–66, feb 2006.  ISSN 1089-778X.  doi: 10.1109/TEVC.2005.851274.  URL http://ieeexplore.ieee.org/document/1583627/.
  * Kohavi and John [2002] R. Kohavi and G. H. John.  Wrappers for feature subset selection.  _Artificial Intelligence_ , 97(1-2):273–324, 2002\.  ISSN 00043702.
  * Li et al. [2013] R. Li, M. T. Emmerich, J. Eggermont, T. Bäck, M. Schütz, J. Dijkstra, and J. H. Reiber.  Mixed integer evolution strategies for parameter optimization.  _Evolutionary Computation_ , 21(1):29–64, 2013\.  ISSN 10636560.
  * Liu et al. [2018] H. Liu, K. Simonyan, and Y. Yang.  Darts: Differentiable architecture search.  _arXiv preprint arXiv:1806.09055_ , 2018.
  * Pfisterer et al. [2019] F. Pfisterer, S. Coors, J. Thomas, and B. Bischl.  Multi-objective automatic machine learning with autoxgboostmc.  _arXiv preprint arXiv:1908.10796_ , 2019.
  * Snoek et al. [2012] J. Snoek, H. Larochelle, and R. P. Adams.  Practical Bayesian Optimization of Machine Learning Algorithms.  In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems 25_ , pages 2951—-2959, 2012.
  * Vanschoren et al. [2014] J. Vanschoren, J. N. van Rijn, B. Bischl, and L. Torgo.  OpenML: networked science in machine learning.  _ACM SIGKDD Explorations Newsletter_ , 15(2):49–60, 2014.  ISSN 19310145.  URL http://dl.acm.org/citation.cfm?doid=2641190.2641198.
  * Waqas et al. [2009] K. Waqas, R. Baig, and S. Ali.  Feature subset selection using multi-objective genetic algorithms.  In _IEEE 13th International Multitopic Conference_ , pages 1–6, 2009\.  ISBN 9781424448722.
  * Xue et al. [2013] B. Xue, M. Zhang, and W. N. Browne.  Particle swarm optimization for feature selection in classification: A multi-objective approach.  _IEEE Transactions on Cybernetics_ , 43(6):1656–1671, 2013.  ISSN 21682267.
  * Xue et al. [2016] B. Xue, M. Zhang, W. N. Browne, and X. Yao.  A Survey on Evolutionary Computation Approaches to Feature Selection.  _IEEE Transactions on Evolutionary Computation_ , 20(4):606–626, 2016.  ISSN 1089778X.
  * Yang [1999] H. Yang.  Data visualization and feature selection: New algorithms for nongaussian data.  _Advances in Neural Information Processing Systems_ , pages 687—-693, 1999.
  * Yu et al. [2002] K. Yu, L. Ji, and X. Zhang.  Kernel Nearest-Neighbor Algorithm.  _Neural Processing Letters_ , 15(2):147–156, 2002\.  ISSN 13704621.  doi: 10.1023/A:1015244902967.  URL http://link.springer.com/10.1023/A:1015244902967.
  * Zitzler and Thiele [1998] E. Zitzler and L. Thiele.  Multiobjective Optimization Using Evolutionary Algorithms — a Comparative Case Ctudy.  In A. E. Eiben, T. Bäck, M. Schoenauer, and H.-P. Schwefel, editors, _Parallel Problem Solving from Nature — PPSN V_ , pages 292–301, Berlin, Heidelberg, 1998. Springer Berlin Heidelberg.  ISBN 978-3-540-49672-4.

## Supplementary Material

Figure 7: Cluster-dendrogram of filter values. We used a range of filter
scoring methods included in the mlr R-package
(https://cran.r-project.org/web/packages/mlr/index.html) and evaluated each of
them twice on a range of datasets. Distances between filters was calculated as
the average (over datasets) L1-distance between filter-values that were rank-
transformed and scaled to $[0,1]$. The vertical axis gives the distance
between clusters. For deterministic filter methods, both evaluation instances
have distance $0$. The two instances of univariate.model.score differ more
from each other than many other groups, so this method was excluded.

### NSGA-II Ablation Study

We show that our methods introduced in Section 4 have beneficial effect on
optimization performance. For this we perform an ablation study of running a
basic NSGA-II, as well as different versions that have incrementally more
specialized operations included. Figure 8 shows the performance of different
methods that successively include uniformly distributed and geometric
initialization (Section 4.3.1), Hamming-weight preserving mutation (Section
4.3.3), and Filter-ensemble based initialization and mutation (Sections 4.3.2
and 4.3.4). The basic algorithm, closely corresponding to the method shown in
Bouraoui et al. [3], has performance far inferior to all GA methods because
its initial population does not cover the $\mathrm{ffrac}$ dimension well
(Figure 2).

Table 5: Description of different GA runs shown in Figure 8. “Naive” Bernoulli initialization consists of drawing each feature selection bit independently from a $\frac{1}{2}$-Bernoulli distribution. Variant | Initialization | Hamming-weight | Filter-ensemble | Filter-ensemble
---|---|---|---|---
| (Section 4.3.1) | preserving mutation | initialization (Section 4.3.2) | mutation
|  | (Section 4.3.3) |  | (Section 4.3.4)
(1) | Bernoulli (naive) | No | No | No
(2) | uniform | No | No | No
(3) | geom | No | No | No
(4) | geom | Yes | No | No
(5) (GA-MO) | geom | Yes | Yes | No
(6) (GA-MO-FE) | geom | Yes | Yes | Yes
Figure 8: Global rank analysis to study the effect of the different
algorithmic components (Table 5) of the GA based on
$\textrm{domHV}_{\mathrm{gen}}$. Ranks are computed per dataset and algorithm
(ties are ranked by their average rank) and then averaged. Higher values are
better.
Figure 9: Comparison of multi-objective and single-objective methods applied
to the kknn learning algorithm. Shown is both the performance
($\mathrm{mmce}_{\mathrm{test}}$) and the fraction of included features
($\mathrm{ffrac}$) of individuals found after 2000 evaluations by the single-
objective baseline BO-SO (tail end of arrows), BO-MO-FE (head of solid arrows)
and GA-MO-FE (head of dashed arrows). Each dataset (Table 2) is shown, values
are averaged over 10 outer cross-validation runs. Individuals of MO methods
were chosen as described in Section 6.2.
Figure 10: Comparison of multi-objective and single-objective methods applied to the xgboost learning algorithm. Shown is both the performance $\mathrm{mmce}_{\mathrm{test}}$ and the fraction of included features $\mathrm{ffrac}$ of individuals found after 2000 evaluations by the single-objective baseline BO-SO (tail end of arrows), BO-MO-FE (head of solid arrows) and GA-MO-FE (head of dashed arrows). Each dataset (Table 2) is shown, values are averaged over 10 outer cross-validation runs. Individuals of MO methods were chosen as described in Section 6.2. Table 6: Methods that perform hyperparameter tuning and feature selection simultaneously (BO-MO-FE, GA-MO-FE) are compared to similar methods that do not perform the tasks simultaneously. The table shows the $\textrm{domHV}_{\mathrm{gen}}$ after $2000$ evaluations. | Learner | Problem | BO-MO-FE | BO-MO-FE-NJ | GA-MO-FE | GA-MO-FE-NJ
---|---|---|---|---|---|---
1 | kknn | AP_Breast_Colon | 0.9761 | 0.9849 | 0.9761 | 0.9761
2 | kknn | AP_Colon_Kidney | 0.9853 | 0.9595 | 0.9836 | 0.9872
3 | kknn | arcene | 0.9030 | 0.8382 | 0.8549 | 0.8898
4 | kknn | cnae-9 | 0.9879 | 0.9692 | 0.9754 | 0.9662
5 | kknn | hill-valley | 0.6653 | 0.6548 | 0.6490 | 0.6722
6 | kknn | ionosphere | 0.9252 | 0.8581 | 0.9143 | 0.8989
7 | kknn | isolet | 0.9957 | 0.9908 | 0.9925 | 0.9959
8 | kknn | lsvt | 0.8856 | 0.8675 | 0.9027 | 0.8770
9 | kknn | madeline | 0.9128 | 0.7505 | 0.9027 | 0.8894
10 | kknn | madelon | 0.9047 | 0.6201 | 0.8881 | 0.8686
11 | kknn | semeion | 0.9872 | 0.9606 | 0.9817 | 0.9805
12 | kknn | sonar | 0.8683 | 0.8374 | 0.8338 | 0.8881
13 | kknn | tecator | 0.9458 | 0.9162 | 0.9399 | 0.9368
14 | kknn | wdbc | 0.9610 | 0.9297 | 0.9565 | 0.9614
15 | SVM | AP_Breast_Colon | 0.9777 | 0.9401 | 0.9698 | 0.9777
16 | SVM | AP_Colon_Kidney | 0.9835 | 0.9595 | 0.9780 | 0.9107
17 | SVM | arcene | 0.9377 | 0.8732 | 0.8898 | 0.7199
18 | SVM | cnae-9 | 0.9862 | 0.9660 | 0.9597 | 0.9590
19 | SVM | hill-valley | 0.7525 | 0.7093 | 0.7599 | 0.7653
20 | SVM | ionosphere | 0.9366 | 0.9082 | 0.9270 | 0.9318
21 | SVM | isolet | 0.9941 | 0.9941 | 0.9942 | 0.9942
22 | SVM | lsvt | 0.8887 | 0.8863 | 0.8861 | 0.8821
23 | SVM | madeline | 0.9205 | 0.6468 | 0.9138 | 0.7789
24 | SVM | madelon | 0.9048 | 0.6219 | 0.8992 | 0.7389
25 | SVM | semeion | 0.9849 | 0.9756 | 0.9814 | 0.9817
26 | SVM | sonar | 0.8718 | 0.8442 | 0.8599 | 0.8847
27 | SVM | tecator | 0.9624 | 0.9469 | 0.9587 | 0.9633
28 | SVM | wdbc | 0.9619 | 0.9339 | 0.9600 | 0.9599
29 | xgboost | AP_Breast_Colon | 0.9777 | 0.9738 | 0.9761 | 0.9777
30 | xgboost | AP_Colon_Kidney | 0.9799 | 0.9800 | 0.9853 | 0.9781
31 | xgboost | arcene | 0.8494 | 0.8405 | 0.8249 | 0.8349
32 | xgboost | cnae-9 | 0.9555 | 0.9434 | 0.9473 | 0.9546
33 | xgboost | hill-valley | 0.7205 | 0.6932 | 0.7510 | 0.7390
34 | xgboost | ionosphere | 0.9203 | 0.8928 | 0.9232 | 0.9239
35 | xgboost | isolet | 0.9909 | 0.9877 | 0.9908 | 0.9892
36 | xgboost | lsvt | 0.8952 | 0.8376 | 0.8945 | 0.8559
37 | xgboost | madeline | 0.8994 | 0.8502 | 0.8915 | 0.8906
38 | xgboost | madelon | 0.8950 | 0.8405 | 0.8921 | 0.8847
39 | xgboost | semeion | 0.9778 | 0.9717 | 0.9817 | 0.9845
40 | xgboost | sonar | 0.8534 | 0.8245 | 0.8615 | 0.8378
41 | xgboost | tecator | 0.9564 | 0.9497 | 0.9583 | 0.9654
42 | xgboost | wdbc | 0.9576 | 0.9310 | 0.9533 | 0.9529