# Optimizing the diffusion coefficient of overdamped Langevin dynamics

T. Lelièvre1, G. Pavliotis2, G. Robin3, R. Santet1 and G. Stoltz1
1: CERMICS, École des Ponts, Marne-la-Vallée, France & MATHERIALS project-
team, Inria Paris, France
2: Department of Mathematics, Imperial College, London, United-Kingdom
3: CNRS & Laboratoire de Mathématiques et Modélisation d’Évry, Évry, France

###### Abstract

Overdamped Langevin dynamics are reversible stochastic differential equations
which are commonly used to sample probability measures in high-dimensional
spaces, such as the ones appearing in computational statistical physics and
Bayesian inference. By varying the diffusion coefficient, there are in fact
infinitely many overdamped Langevin dynamics which are reversible with respect
to the target probability measure at hand. This suggests to optimize the
diffusion coefficient in order to increase the convergence rate of the
dynamics, as measured by the spectral gap of the generator associated with the
stochastic differential equation. We analytically study this problem here,
obtaining in particular necessary conditions on the optimal diffusion
coefficient. We also derive an explicit expression of the optimal diffusion in
some appropriate homogenized limit. Numerical results, both relying on
discretizations of the spectral gap problem and Monte Carlo simulations of the
stochastic dynamics, demonstrate the increased quality of the sampling arising
from an appropriate choice of the diffusion coefficient.

##  1 Introduction

Predicting properties of materials and macroscopic physical systems in the
framework of statistical physics [8], using for instance molecular dynamics
[30, 68, 44, 49, 5] and obtaining the distribution of parameter values in
Bayesian inference [60], both require the sampling of probability measures in
high-dimensional spaces. Methods of choice to sample such probability measures
rely on stochastic dynamics, in particular Markov Chain Monte Carlo (MCMC)
methods. The convergence of these methods may however be quite slow because
the target measure is typically concentrated on a few high probability modes
separated by low probability regions, or because the Hessian of the
corresponding potential, see (1), is badly conditioned. Various algorithms
have been proposed over the years to improve the sampling of such probability
measures, for instance using importance sampling strategies, and/or
interacting replicas; see for instance the extensive review [41] in the
context of molecular dynamics.

##### Sampling Boltzmann–Gibbs measures.

We focus in this work on overdamped Langevin dynamics, a popular choice in
molecular dynamics, computational statistics [62], and also machine learning
(see for instance [20]). This dynamics is ergodic (and, in fact, reversible)
with respect to the Boltzmann–Gibbs distribution, which is the probability
measure with density

$\mu(q)=Z^{-1}\mathrm{e}^{-V(q)},\qquad
Z=\int_{\mathcal{Q}}\mathrm{e}^{-V}<+\infty,$ (1)

where $V\in\mathcal{C}^{\infty}(\mathcal{Q})$. In molecular dynamics, the
measure $\mu$ is the distribution sampled by the possible states of a system
in the canonical thermodynamic ensemble. In that case, the system is at
thermal equilibrium with a heat bath of fixed temperature $T$, which is taken
into account by considering $\beta V$ instead of $V$ in (1) with
$\beta=1/(k_{\rm B}T)$, where $k_{\rm B}$ is the Boltzmann constant (in this
work, we simply set $\beta=1$). To simplify the presentation, and because
periodic boundary conditions are commonly used for molecular dynamics
simulations, we restrict ourselves to the case when the configuration space
$\mathcal{Q}$ is the $d$-dimensional torus $\mathbb{T}^{d}$ (with
$\mathbb{T}=\mathbb{R}\backslash\mathbb{Z}$ the one dimensional torus). Many
arguments in our work can be extended, at least formally, to systems defined
in unbounded configuration spaces.

The overdamped Langevin dynamics is given by the following stochastic
differential equation (SDE):

$dq_{t}=-\nabla V(q_{t})\,dt+\sqrt{2}\,dW_{t},$ (2)

where $(W_{t})_{t\geqslant 0}$ is a standard $d$-dimensional Brownian motion.
This dynamics enjoys many nice properties. In particular, it is reversible
with respect to the Gibbs measure (1), and, under appropriate s on the
potential energy function $V$ (see for instance [7]), its marginal law in time
converges exponentially fast to the target distribution (1) in various norms
and distances, e.g. total variation, relative entropy (Kullback–Leibler),
$L^{2}(\mu)$, or weighted $L^{\infty}$ norms. For example, when the overdamped
Langevin dynamics in $\mathbb{R}^{d}$ is considered, it is sufficient for the
potential $V$ to be smooth and convex at infinity to ensure the existence of
Poincaré and logarithmic Sobolev inequalities for the Gibbs measure which, in
turn, leads to exponentially fast convergence to the target distribution in
$L^{2}(\mu)$ and in relative entropy, respectively. This class of potentials
includes many potentials of interest, for example multistable potentials.
Naturally, in the nonconvex case, the convergence rate typically scales very
badly with respect to the dimension $d$ (and the temperature). In this work,
we will consider convergence in $L^{2}(\mu)$, _i.e._ the
$\chi^{2}$-divergence. This is the relevant setting _e.g._ for obtaining error
estimates on trajectory averages, relying on the Central Limit Theorem.

It is important to note that the constant in the Poincaré and in the
logarithmic Sobolev inequalities is a property of the Gibbs measure. In
particular, the rate of convergence to the target measure for the overdamped
Langevin dynamics (2) depends only on the potential $V$, and convergence can
be very slow when entropic or energetic barriers are present. It is therefore
natural to consider alternatives to the overdamped Langevin dynamics that
include “hyperparameters” that can be tuned to improve the rate of convergence
to equilibrium and reduce the asymptotic variance.

There are in fact infinitely many overdamped Langevin dynamics which admit (1)
as invariant probability measure; see [54, Section 4.6]. A comprehensive list
of stochastic dynamics that can be used in order to sample from a Gibbs
measure can be found in [25]. For instance, dynamics of the form

$dq_{t}=\left(-\nabla V(q_{t})+\gamma(q_{t})\right)dt+\sqrt{2}\,dW_{t},$

where the vector field $\gamma:\mathbb{T}^{d}\to\mathbb{R}^{d}$ is smooth and
such that

$\operatorname{div}\left(\gamma\mu\right)=0,$

leave invariant the target probability measure with density (1); and is
ergodic, but not reversible, with respect to this probability measure under
certain s on $V,\gamma$. The hope is that the added drift term $\gamma$ that
renders the dynamics nonreversible, can accelerate convergence, measured in
the $L^{2}(\mu)$ sense, and also reduce the asymptotic variance. This approach
was extensively studied in the past decade [39, 40, 46, 58, 59, 24, 26]. For
example, for Gaussian targets, the optimal nonreversible perturbation $\gamma$
can be obtained in an algorithmic manner [46]. The presence of the divergence-
drift term $\gamma$ can lead to computational difficulties, since the
resulting dynamics can be stiff. Appropriate numerical methods need to be
developed to address this issue [24, 26]. In fact, the discretized
nonreversible dynamics can perform worse than the reversible dynamics if the
stiffness issue is not addressed adequately [71, 51].

##### Optimizing reversible dynamics for better sampling.

Another alternative, which is the focus of this work, is to stay within the
class of reversible dynamics, but with a space-dependent diffusion matrix. For
a given diffusion matrix
$\mathcal{D}\in\mathcal{C}^{1}(\mathbb{T}^{d},\mathcal{S}_{d}^{+})$ (with
$\mathcal{S}_{d}^{+}$ the set of real, symmetric, positive matrices of size
$d\times d$), the associated overdamped Langevin dynamics reads

$dq_{t}=\left(-\mathcal{D}(q_{t})\nabla
V(q_{t})+\operatorname{div}\mathcal{D}(q_{t})\right)dt+\sqrt{2}\,\mathcal{D}^{1/2}(q_{t})\,dW_{t},$
(3)

where $\mathcal{D}^{1/2}$ is defined by functional calculus, and
$\operatorname{div}\,\mathcal{D}$ is the vector whose $i$-th component is the
divergence of the $i$-th column of the matrix
$\mathcal{D}=[\mathcal{D}_{1},\dots,\mathcal{D}_{d}]$, _i.e._
$\operatorname{div}\,\mathcal{D}=\left(\operatorname{div}\,\mathcal{D}_{1},\dots,\operatorname{div}\,\mathcal{D}_{d}\right)^{\top}$.
Modulating the diffusion to improve the sampling was explored to some extent
both in the computational statistics literature [61] and in molecular dynamics
[59, 1]; and also for simulated annealing [28].

The goal of this paper is to optimize the diffusion matrix in order to
accelerate convergence. Intuitively, it seems relevant to accelerate the
diffusion in regions of low probability under the target measure
(equivalently, regions of high values of the energy function $V$), in order to
more efficiently and rapidly find another mode, and slow down the diffusion in
regions of high probability since these are the zones where sampling should be
favored. Clearly, in order for this optimization problem to make sense,
appropriate constraints on the magnitude of the diffusion matrix need to be
introduced. These constraints on the diffusion matrix and their effects on the
solution of the associated constrained optimization problem are discussed in
detail later on.

The $L^{2}(\mu)$ convergence rate of the overdamped Langevin diffusion (3) is
given by the spectral gap of the infinitesimal generator. The latter operator
acts on test functions $\varphi$ as

$\mathcal{L}_{\mathcal{D}}\varphi=\left(-\mathcal{D}\nabla
V+\operatorname{div}(\mathcal{D})\right)^{\top}\nabla\varphi+\mathcal{D}:\nabla^{2}\varphi.$
(4)

It is standard to show (see, _e.g._ [46, 7, 49]) that,
$\mathcal{L}_{\mathcal{D}}$ admits a positive spectral gap
$\Lambda(\mathcal{D})$ if and only if, for any initial condition
$f(0)=\mu_{0}/\mu\in L^{2}(\mu)$, and denoting by $f(t)\mu$ the probability
density function of the process $q_{t}$ at time $t$,

$\forall t\geqslant 0,\,\left\lVert
f(t)-1\right\rVert_{L^{2}(\mu)}\leqslant\mathrm{e}^{-\Lambda(\mathcal{D})t}\left\lVert
f(0)-1\right\rVert_{L^{2}(\mu)},$ (5)

where $\left\lVert\cdot\right\rVert_{L^{2}(\mu)}$ is the norm on $L^{2}(\mu)$,
see Section 2.1 for precise definitions. For a fixed target probability
density $\mu$ (_i.e._ a fixed potential energy function $V$), the generator
$\mathcal{L}_{\mathcal{D}}$, and thus the spectral gap $\Lambda(\mathcal{D})$,
are parameterized by the diffusion matrix $\mathcal{D}$.

The aim of this work is to compute, explicitly or numerically, the optimal
diffusion matrix leading to the largest spectral gap, and thus to the largest
convergence rate. For reversible Markov chains on discrete spaces, this
question was explored in [16, 15]. A subtle issue in this endeavour is the
normalization of the diffusion matrix $\mathcal{D}$. Indeed, the convergence
rate in (5) can trivially be increased by a factor $\alpha\geqslant 1$ upon
multiplying $\mathcal{D}$ by $\alpha$. In fact, it seems advantageous to let
$\alpha$ go to infinity. The catch is of course that one should compare
dynamics which are defined on similar timescales. From a practical point of
view, one way of setting a timescale is to consider discretizations of
overdamped Langevin dynamics. Larger diffusion matrices require smaller values
of the time step. The issue of normalization is discussed more thoroughly in
Section 2.2.

The optimization of the diffusion matrix is related to a body of literature on
accelerating the convergence of Langevin-like dynamics through particular
choices for the diffusion matrix $\mathcal{D}$. This is very much related to
the optimization of the damping coefficient [19] or of the mass matrix in the
Langevin Monte Carlo method, a.k.a. as preconditioning, see [18, 10, 32, 45,
34, 21]. One example of such methods is the Riemannian manifold Langevin Monte
Carlo method, introduced in [32], and which reduces to (3) in the so-called
overdamped limit. There, for strictly convex potentials, the diffusion matrix
$\mathcal{D}$ is chosen as the inverse of the Fisher–Rao metric tensor
associated with the target measure, namely the inverse of the Hessian of $V$.
From a functional analysis point of view, this can be understood, for a
problem set on $\mathbb{R}^{d}$, from the Brascamb-Lieb inequality [37, 64]:

$\mbox{Var}_{\mu}(h)\leqslant\int_{\mathbb{R}^{d}}\langle\nabla
h,(\mbox{Hess}(V))^{-1}\nabla h\rangle\,\mu(q)\,dq,$

for an arbitrary observable $h\in\mathcal{C}^{1}$ with
$\mbox{Var}_{\mu}(h)<+\infty$. In particular, when the diffusion matrix is
chosen to be the inverse of the Hessian of $V$, this leads to a spectral gap
for the dynamics (3) that is independent of $V$ and of the number of
dimensions and is thus “optimal” in this sense. We emphasize that, although
relevant, this unconstrained optimization problem for the diffusion matrix is
quite different from the one we consider in this paper: we optimize the
spectral gap for the dynamics on the torus, with no convexity s on $V$ and
with constraints on the magnitude of $\mathcal{D}$.

At the time we completed this work, we became aware of the preprint [23] where
the authors consider a very similar problem, that they solve using different
techniques, namely Stein Kernels and moment measures. More precisely, by using
an appropriate transportation map, they transform the nonconvex log-
probability to a convex one, for which the Brascamp–Lieb inequality yields the
optimal diffusion, for a well-chosen normalizing constraint. We intend to
further explore the links between this recent preprint and ours in future
works.

Let us also mention results obtained via the theory of partial differential
equations for eigenvalue problems associated with operators of the form
$-\operatorname{div}(\mathcal{D}\nabla\cdot)$, see [38, Chapter 10]. The main
difference with such results and ours is that we consider the Laplace
operators on the weighted space $L^{2}(\mu)$ as we are interested in sampling
the nontrivial probability measure $\mu$.

##### Main contributions and extensions.

The main contributions of this work are the following.

  * •

First, we formulate the optimization of the convergence rate of the Langevin
dynamics (3) with respect to $\mathcal{D}$ as a convex optimization program,
for which the well-posedness is guaranteed by a simple adaptation of standard
results in the literature [38]. We also discuss some important properties of
the resulting optimal diffusion matrix, such as its formal characterization
using the Euler–Lagrange equation as well as its positivity. Note that, for
Markov chains, it was already observed in [16] that maximizing the spectral
gap can be formulated as a convex optimization problem.

  * •

Second, we propose, in low dimensional scenarios, a numerical procedure
combining a finite element parametrization and an optimization algorithm, to
compute the optimal diffusion matrix in practice. The numerical procedure is
illustrated on simple one-dimensional examples. All the methods and
experiments are provided in open source Python and Julia codes, available at
https://github.com/rsantet/Optimal_Overdamped_Langevin_Diffusion_Python and
https://github.com/rsantet/Optimal_Overdamped_Langevin_Diffusion_Julia.

  * •

Our third main contribution is to study the behavior of the optimal diffusion
matrix in the homogenized limit where the target measure is associated with a
highly oscillating potential. In a one-dimensional setting, we show that the
optimal diffusion matrix has an analytical expression, proportional to the
inverse of the target density, which is in accordance with various previous
heuristics (as in [61]). This analytical solution can be used as an
initialization guess for the optimization algorithm mentioned above, or as a
proxy for the optimal diffusion matrix which does not require costly convex
optimization procedures.

  * •

Our last main contribution is to propose a sampling algorithm based on a
simple Random Walk Metropolis Hastings algorithm with a proposal variance
depending on the current state, which takes advantage of the precomputed
optimal diffusion matrix. We show that the associated Markov chain converges
to the Langevin dynamics with optimal diffusion matrix in the diffusive limit.
We also present some numerical experiments illustrating the behavior of the
complete procedure on simple examples, and highlighting in particular that the
dynamics is less metastable when using the optimal diffusion matrix or its
homogenized approximation. We show that the homogenization approximation
yields similar results as the optimal diffusion matrix, which makes it a good
choice in practice since it has an analytical expression.

We believe that this work calls for various perspectives and extensions, in
particular:

  * •

A key issue remains to find a relevant normalization constraint e.g. based on
a numerical criterion, as discussed in Remark 3;

  * •

In various contexts, we restrict ourselves to isotropic diffusion matrices or
even to the one-dimensional setting (in particular for numerical simulations).
Considering a larger class of positive definite symmetric matrices may be
important to tackle anisotropic targets;

  * •

In order to apply the methodology to actual systems of interest, one needs to
adapt our numerical approach to higher dimensional scenarios, even in the case
of a scalar diffusion coefficient, i.e. when the diffusion matrix is a
multiple of the identity. It is indeed impossible to numerically optimize a
diffusion matrix which genuinely depends on the full configuration $q$ using
discretization techniques such as finite elements, because of the curse of
dimensionality of the corresponding approximation space. Apart from resorting
to parsimonious functional representations of the diffusion matrix (which
usually presuppose some form of smoothness), it seems appropriate to look for
diffusion matrices parameterized by a low dimensional function $\xi(q)$, a so-
called reaction coordinate in molecular dynamics [47], which summarizes some
important information on the system. In this context, the potential energy
function is replaced by the free energy associated with $\xi$, and the
optimization problem can again be reformulated in a low dimensional setting
where finite element methods can be used.

  * •

In this paper we consider the problem of maximizing the spectral gap. The
self-adjointness of the generator implies that maximizing the spectral gap is
typically equivalent to minimizing the asymptotic variance, uniformly over the
space of $L^{2}(\mu)$ observables [6]. We will explore the combination of
spectral gap maximization and asymptotic variance minimization in a future
work. In particular, our aim will be to minimize an appropriate cost
functional that measures the true cost of the actual computation, by using the
approach developed in [33].

  * •

The objective in this paper is to optimize the convergence in
$\chi^{2}$-divergence (see (5)). A natural question would be to consider
convergence in other distances, such as relative entropy with respect to the
target distribution (which amounts to optimizing Logarithmic Sobolev
Inequality constants).

##### Outline of the work.

We formally define the maximization of the spectral gap in Section 2 and make
precise the normalization of the diffusion matrix. We then study the well
posedness of the maximization problem in Section 3, and characterize the
optimal diffusion matrix. We next present in Section 4 how to numerically
solve the problem, using in particular a finite element method to discretize
the diffusion matrix. As an alternative to a full-scale numerical simulations,
or in order to start the optimization procedure with a relevant initial guess,
we consider in Section 5 the optimal diffusion matrix arising in the
homogenization limit where the given probability density is periodically
replicated with a decreasing spatial period. We finally demonstrate in Section
6 with Monte Carlo simulations that an optimized diffusion matrix is
beneficial for the convergence of the dynamics. In the appendices we present
the proofs of various technical results.

##  2 Formulation of the optimization problem

We present in this section the problem of interest, namely the maximization of
the spectral gap of the generator of the Langevin dynamics (3) with invariant
measure $\mu$ given by (1). We start by precisely formulating the optimization
problem in Section 2.1. Then, in Section 2.2, we discuss normalization
constraints which are required for the well-posedness of the problem. For
simplicity, we restrict ourselves to the case when the diffusion process under
consideration takes as input a configuration in the $d$-dimensional torus
$\mathbb{T}^{d}$. Various results can easily be adapted to more general
domains.

###  2.1 Maximizing the spectral gap of Langevin dynamics

As discussed in the introduction around (5), the rate of convergence of the
dynamics (3) towards the invariant measure $\mu$ is governed by the spectral
gap of the operator $\mathcal{L}_{\mathcal{D}}$ (see _e.g._ [7, Chapter 4]).
In order to precisely define this quantity in our context, we work on the
Hilbert space $L^{2}(\mu)$. For a given diffusion matrix
$\mathcal{D}:\mathbb{T}^{d}\to\mathcal{S}_{d}^{+}$, a simple computation shows
that the generator $\mathcal{L}_{\mathcal{D}}$, considered as an operator on
$L^{2}(\mu)$, can be written as

$\mathcal{L}_{\mathcal{D}}=-\nabla^{*}\mathcal{D}\nabla=-\sum_{i,j=1}^{d}\partial_{q_{j}}^{*}\mathcal{D}_{i,j}\partial_{q_{i}},$
(6)

where $A^{*}$ is the $L^{2}(\mu)$-adjoint of a closed operator $A$. In
particular, $\partial_{q_{i}}^{*}=-\partial_{q_{i}}+\partial_{q_{i}}V$. The
quadratic form associated with $\mathcal{L}_{\mathcal{D}}$ is easily obtained
from the expression (6) of the generator: formally, for any
$u:\mathbb{T}^{d}\to\mathbb{R}$,

$\left\langle
u,\mathcal{L}_{\mathcal{D}}u\right\rangle_{L^{2}(\mu)}=-\int_{\mathbb{T}^{d}}\nabla
u(q)^{\top}\mathcal{D}(q)\nabla u(q)\,\mu(q)\,dq.$

The right-hand side of the above equality is nonpositive, but possibly
infinite at this level of generality. It is finite for instance when
$\mathcal{D}\in L^{\infty}(\mathbb{T}^{d},\mathcal{S}_{d}^{+})$ and $u\in
H^{1}(\mu)$, the subspace of $L^{2}(\mu)$ composed of functions in
$L^{2}(\mu)$ whose (distributional) partial derivatives with respect to
$q_{i}$ also belong to $L^{2}(\mu)$. We restrict for simplicity to this
setting in the sequel, and consider in particular diffusion matrices with
bounded coefficients.

Since 0 is an eigenvalue of $\mathcal{L}_{\mathcal{D}}$ whose associated
eigenvectors are constant functions, we introduce the following subspace to
define the spectral gap:

$H^{1}_{0}(\mu)=\left\\{u\in H^{1}(\mu)\ \middle|\
\int_{\mathbb{T}^{d}}u(q)\,\mu(q)\,dq=0\right\\}.$

We also define $H^{1}_{0}(\mathbb{T}^{d})$ as the space
$H^{1}_{0}(\mathbf{1})$. The spectral gap is finally defined as

$\Lambda(\mathcal{D})=\inf_{u\in
H^{1}_{0}(\mu)\setminus\\{0\\}}\frac{\displaystyle\int_{\mathbb{T}^{d}}\nabla
u(q)^{\top}\mathcal{D}(q)\nabla
u(q)\,\mu(q)\,dq}{\displaystyle\int_{\mathbb{T}^{d}}u(q)^{2}\mu(q)\,dq}.$ (7)

For a bounded positive definite diffusion matrix $\mathcal{D}$ satisfying
$\mathcal{D}(q)\geqslant c\mathrm{Id}_{d}$ almost everywhere for some $c>0$
(the inequality being understood in the sense of symmetric matrices), standard
arguments from calculus of variations show that the infimum is attained for an
eigenvector $u_{\mathcal{D}}$ associated with the second eigenvalue of
$-\mathcal{L}_{\mathcal{D}}$ (ranked by increasing values), _i.e._ with the
first nonzero eigenvalue of $-\mathcal{L}_{\mathcal{D}}$; and that
$\Lambda(\mathcal{D})$ is in fact the first nonzero eigenvalue of
$-\mathcal{L}_{\mathcal{D}}$ (see for instance [38, Chapter 1]).

######  Remark 1 (Degenerate diffusion matrices).

If there exists $q\in\mathbb{T}^{d}$ such that $\mathcal{D}(q)$ is rank
deficient (_i.e._ $\xi^{\top}\mathcal{D}(q)\xi=0$ for some
$\xi\in\mathbb{R}^{d}\setminus\\{0\\}$), the domain of the quadratic form
associated with $\mathcal{L}_{\mathcal{D}}$ can in fact be larger than
$H^{1}_{0}(\mu)$; see [72]. In this case, $\Lambda(\mathcal{D})$ is in general
no longer equal to the spectral gap of the operator
$\mathcal{L}_{\mathcal{D}}$, but is larger than the spectral gap (and hence
provides a too optimistic convergence rate) since the infimum is taken on a
set of functions not sufficiently large. See [38, Remark 1.2.3] for a
discussion on sufficient conditions for the operator
$\mathcal{L}_{\mathcal{D}}$ to have a discrete spectrum even in the presence
of degeneracies of the diffusion matrix.

######  Remark 2 (Unbounded diffusion matrices).

We emphasize the fact that we will restrict ourselves to bounded diffusion
matrices; otherwise the upper semicontinuity of the mapping
$\mathcal{D}\mapsto\Lambda(\mathcal{D})$ is not clear (see Lemma 3). The
latter property is crucial to ensure that the maximization problem (8) below
is well posed.

We propose to maximize the spectral gap of the overdamped Langevin dynamics
with respect to the diffusion matrix $\mathcal{D}$. This corresponds here to
solving the following optimization problem:

$\Lambda^{\star}=\sup_{\mathcal{D}\in\mathfrak{D}}\Lambda(\mathcal{D}),$ (8)

where $\mathfrak{D}$ is a subset of bounded measurable functions from
$\mathbb{T}^{d}$ to $\mathcal{S}_{d}^{+}$. The choices we consider in the
sequel for $\mathfrak{D}$ take into account that the diffusion needs to be
normalized in some way since, for any
$\mathcal{D}:\mathbb{T}^{d}\rightarrow\mathcal{S}_{d}^{+}$,

$\forall t\geqslant 0,\qquad\Lambda(t\mathcal{D})=t\Lambda(\mathcal{D}),$ (9)

so that a maximization over a set containing a line $\mathbb{R}\mathcal{D}$
with $\Lambda(\mathcal{D})>0$ would lead to $\Lambda^{\star}=+\infty$.

######  Remark 3 (A numerical motivation for normalizing the diffusion
matrix).

The need for a normalization can also be motivated by numerical
considerations. For instance, an Euler–Maruyama discretization of (3) with a
time step $\Delta t>0$ reads

$q^{n+1}=q^{n}+\Delta t\left[-\mathcal{D}\nabla
V+\operatorname{div}\,\mathcal{D}\right](q^{n})+\sqrt{2\Delta
t}\,\mathcal{D}^{1/2}(q^{n})\,G^{n},$ (10)

where $(G^{n})_{n\geqslant 0}$ is a family of independent and identically
distributed standard $d$-dimensional Gaussian random variables. It is apparent
from the above formula that the behavior of the numerical scheme actually
depends on $\Delta t\,\mathcal{D}$ (and the divergence of this quantity), and
not on $\mathcal{D}$ or $\Delta t$ alone, so that multiplying $\mathcal{D}$ by
a constant can be equivalently seen as multiplying the time step by the same
constant. In this context, one should think of normalization criteria as ways
to fix the amount of numerical error involved in the discretization, for
instance by fixing the average rejection probability in a Metropolis–Hastings
scheme whose proposal is provided by (10). The normalization we discuss in
Section 2.2 is conceptually simpler than such a criterion, and not related to
a specific numerical scheme. We intend to further study this perspective in
following works.

###  2.2 Normalization constraint

We make precise in this section the normalization we consider for diffusion
matrices, in view of the discussion after (8). In essence, we choose a
(weighted) Lebesgue norm to be smaller than 1.

A quantity which naturally appears in the mathematical formulation of the
problem is the product $\mathcal{D}(q)\mathrm{e}^{-V(q)}$, present for
instance in the numerator of the Rayleigh ratio (7). Note also that the
divergence of this product is the drift in (3) (up to the factor
$\mathrm{e}^{-V(q)}$). This motivates normalizing the product
$\mathcal{D}\mathrm{e}^{-V}$ (which is proportional to $\mathcal{D}\mu$)
rather than $\mathcal{D}$ itself.

To formulate precisely the normalization we consider, it is convenient to
introduce, for $a,b\geqslant 0$, the set

$\mathcal{M}_{a,b}=\left\\{M\in\mathcal{S}_{d}^{+}\,\middle|\,\forall\xi\in\mathbb{R}^{d},\
a|\xi|^{2}\leqslant\xi^{\top}M\xi\leqslant\frac{1}{b}|\xi|^{2}\right\\},$

with $|\cdot|$ the Euclidean norm on $\mathbb{R}^{d}$, and the convention that
the right-hand side of the last inequality is $+\infty$ when $b=0$.
Considering $a,b>0$ is useful to study homogenization limits, as we do in
Section 5. We will always require that

$ab\leqslant 1,$

in order for $\mathcal{M}_{a,b}$ to be nonempty. We endow $\mathcal{M}_{a,b}$
with a norm $|\cdot|_{\rm F}$, which is any matrix norm compatible with the
order on symmetric positive matrices, _i.e._
$\left|M_{1}\right|_{\mathrm{F}}\leqslant\left|M_{2}\right|_{\mathrm{F}}$ when
$0\leqslant M_{1}\leqslant M_{2}$ in the sense of symmetric matrices. One such
norm is the Frobenius norm.

We next denote by $L^{p}_{V}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ for
$p\in[1,+\infty]$ the Banach space of measurable functions $\mathcal{D}$ from
$\mathbb{T}^{d}$ to $\mathcal{S}_{d}^{+}$ such that
$\mathcal{D}(q)\mathrm{e}^{-V(q)}\in\mathcal{M}_{a,b}$ for almost every
$q\in\mathbb{T}^{d}$ and
$\left\lVert\mathcal{D}\right\rVert_{L^{p}_{V}}<+\infty$, with the norm

$\left\lVert\mathcal{D}\right\rVert_{L^{p}_{V}}=\left(\int_{\mathbb{T}^{d}}\left|\mathcal{D}(q)\right|_{\mathrm{F}}^{p}\,\mathrm{e}^{-pV(q)}\,dq\right)^{1/p},$
(11)

when $p<+\infty$, and
$\left\lVert\mathcal{D}\right\rVert_{L^{\infty}_{V}}=\left\|\left|\mathcal{D}\right|_{\mathrm{F}}\mathrm{e}^{-V}\right\|_{L^{\infty}(\mathbb{T}^{d})}$
when $p=+\infty$. Let us emphasize that the $L^{p}_{V}$-norm is in fact the
usual $L^{p}$-norm applied to the product
$\left|\mathcal{D}\right|_{\mathrm{F}}\mathrm{e}^{-V}$. Other choices of
normalization could be considered, see Remark 8 in Section D.4. One motivation
for (11) is that it leads to simple expressions for optimal diffusion matrices
in the homogenized limit (see Section 5.2.2).

Simply considering a normalization constraint in $L^{\infty}_{V}$-norm leads
to a trivial solution for the maximization of the spectral gap (which simply
saturates the constraint), as shown in the following result (see related
discussions in [38, Section 10.1]).

######  Lemma 1.

Fix $a\geqslant 0$ and $b>0$ such that $ab\leqslant 1$. Then, the spectral gap
$\Lambda(\mathcal{D})$ is maximized on
$L^{\infty}_{V}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ by the diffusion matrix
$\mathcal{D}^{\star}(q)=b^{-1}\mathrm{e}^{V}\mathrm{Id}_{d}$.

###### Proof.

Note first that, for any $\mathcal{D}\in
L^{\infty}_{V}(\mathbb{T}^{d},\mathcal{M}_{a,b})$, it holds
$\mathcal{D}(q)\leqslant b^{-1}\mathrm{e}^{V(q)}\mathrm{Id}_{d}$ in the sense
of symmetric matrices for almost all $q\in\mathbb{T}^{d}$. Therefore, for any
$u\in H^{1}_{0}(\mu)$,

$\frac{\displaystyle\int_{\mathbb{T}^{d}}\nabla
u(q)^{\top}\mathcal{D}(q)\nabla
u(q)\,\mu(q)\,dq}{\displaystyle\int_{\mathbb{T}^{d}}u(q)^{2}\,\mu(q)\,dq}\leqslant\frac{\displaystyle\int_{\mathbb{T}^{d}}\nabla
u(q)^{\top}\mathcal{D}^{\star}(q)\nabla
u(q)\,\mu(q)\,\,dq}{\displaystyle\int_{\mathbb{T}^{d}}u(q)^{2}\,\mu(q)\,dq}.$

Passing to the infimum over $u$ on the left-hand side first, then on the
right-hand side, leads to
$\Lambda(\mathcal{D})\leqslant\Lambda(\mathcal{D}^{\star})$ for all
$\mathcal{D}\in L^{\infty}_{V}(\mathbb{T}^{d},\mathcal{M}_{a,b})$. Moreover,
it obviously holds that $b^{-1}\mathrm{e}^{V}\mathrm{Id}_{d}\in
L^{\infty}_{V}(\mathbb{T}^{d},\mathcal{M}_{a,b})$, which concludes the proof.
∎

In view of Lemma 1, we focus in the remainder of this work on a constraint on
the $L^{p}_{V}$ norm of the diffusion matrix for $1\leqslant p<+\infty$, and
therefore consider the subsets

$\mathfrak{D}_{p}^{a,b}=\left\\{\mathcal{D}\in
L^{\infty}_{V}(\mathbb{T}^{d},\mathcal{M}_{a,b})\,\middle|\,\int_{\mathbb{T}^{d}}\left|\mathcal{D}(q)\right|_{\mathrm{F}}^{p}\,\mathrm{e}^{-pV(q)}\,dq\leqslant
1\right\\}.$ (12)

In order for these sets to be nonempty, the parameter $a\geqslant 0$ should be
small enough.

######  Lemma 2.

Fix $1\leqslant p<+\infty$, and $a,b\geqslant 0$. The set
$\mathfrak{D}_{p}^{a,b}$ is nonempty if and only if

$ab\leqslant 1,\qquad a\leqslant
a_{\mathrm{max}}:=\frac{1}{\left|\mathrm{Id}_{d}\right|_{\mathrm{F}}}.$ (13)

###### Proof.

Note that any element of $\mathfrak{D}_{p}^{a,b}$ is lower bounded by
$a\mathrm{e}^{V}\mathrm{Id}_{d}$. Therefore, the set is nonempty if and only
if it contains $a\mathrm{e}^{V}\mathrm{Id}_{d}$. The equality

$\int_{\mathbb{T}^{d}}\left\lvert
a\mathrm{e}^{V}\mathrm{Id}_{d}\right\rvert_{\rm
F}^{p}\mathrm{e}^{-pV}=a^{p}\left\lvert\mathrm{Id}_{d}\right\rvert_{\rm
F}^{p}$

then gives the claimed result. ∎

We can now state the maximization problem which is considered in the remainder
of this work. Fix $p\in[1,+\infty)$. For any $a\in[0,a_{\mathrm{max}}]$ and
$b>0$ such that $ab\leqslant 1$,

$\textrm{find }\mathcal{D}^{\star}\in\mathfrak{D}_{p}^{a,b}\textrm{ such that
}\Lambda(\mathcal{D}^{\star})=\sup_{\mathcal{D}\in\mathfrak{D}_{p}^{a,b}}\Lambda(\mathcal{D}).$
(14)

##  3 Theoretical analysis of the spectral gap optimization

We present in this section some theoretical results on the mathematical
analysis of the optimization problem (14). We start by well-posedness results
in Section 3.1, and we then formally characterize in Section 3.2 the optimal
diffusion matrix.

###  3.1 Well posedness of the optimization problem

The aim of this section is to show that the maximization of the spectral gap
$\Lambda(\mathcal{D})$ on the set $\mathfrak{D}_{p}^{a,b}$ defined in (12) is
a well-posed problem. This is easily done by adapting arguments from [38].
Some results rely on the fact that the probability measure $\mu$ satisfies a
Poincaré inequality. Recall that the potential $V$ is assumed to be
$\mathcal{C}^{\infty}(\mathbb{T}^{d})$ in all this work. This implies in
particular that the potential $V$ is bounded and that the density $\mu$ is
bounded from below by a positive constant and bounded from above on
$\mathbb{T}^{d}$. The measure with density (1) therefore satisfies a Poincaré
inequality, _i.e._ there exists $C_{\mu}>0$ such that

$\inf_{u\in
H^{1}_{0}(\mu)}\left\\{\frac{\displaystyle\int_{\mathbb{T}^{d}}|\nabla
u(q)|^{2}\,\mu(q)\,dq}{\displaystyle\int_{\mathbb{T}^{d}}u(q)^{2}\,\mu(q)\,dq}\right\\}=\frac{1}{C_{\mu}}.$
(15)

In addition, there exists $u_{\mu}\in H^{1}_{0}(\mu)\setminus\\{0\\}$ such
that the infimum (15) is attained, _i.e._

$\frac{\displaystyle\int_{\mathbb{T}^{d}}|\nabla
u_{\mu}(q)|^{2}\,\mu(q)\,dq}{\displaystyle\int_{\mathbb{T}^{d}}u_{\mu}(q)^{2}\,\mu(q)\,dq}=\frac{1}{C_{\mu}}.$

Let us first prove that the function $\Lambda$ is well-defined on
$\mathfrak{D}_{p}^{a,b}$ for $a,b\geqslant 0$. On the one hand,
$\Lambda(\mathcal{D})\geqslant 0$ (since $\mathcal{D}$ is positive semi-
definite on $\mathfrak{D}_{p}^{a,b}$). On the other hand, for a fixed function
$u\in\mathcal{C}^{\infty}(\mathbb{T}^{d})$ satisfying
$\int_{\mathbb{T}^{d}}u(q)\,\mu(q)\,dq=0$ and
$\int_{\mathbb{T}^{d}}u(q)^{2}\mu(q)\,dq=1$, the function $|\nabla u(q)|^{2}$
is bounded on $\mathbb{T}^{d}$. Since there exists a constant
$K\in\mathbb{R}_{+}$ such that $0\leqslant\xi^{\top}\mathcal{D}(q)\xi\leqslant
K\left|\mathcal{D}(q)\right|_{\mathrm{F}}|\xi|^{2}$ (the constant $K$ being 1
for the Frobenius norm), it holds (recall that $Z$ is defined in (1))

$\displaystyle\Lambda(\mathcal{D})$ $\displaystyle\leqslant
Z^{-1}K\left\lVert|\nabla
u|^{2}\right\rVert_{L^{\infty}(\mathbb{T}^{d})}\int_{\mathbb{T}^{d}}\left|\mathcal{D}(q)\right|_{\mathrm{F}}\,\mathrm{e}^{-V(q)}\,dq$
$\displaystyle\leqslant Z^{-1}K\left\lVert|\nabla
u|^{2}\right\rVert_{L^{\infty}(\mathbb{T}^{d})}\left(\int_{\mathbb{T}^{d}}\left|\mathcal{D}(q)\right|_{\mathrm{F}}^{p}\mathrm{e}^{-pV(q)}\,dq\right)^{1/p}<+\infty,$

using the Hölder’s inequality on $\mathbb{T}^{d}$. This shows that
$\Lambda(\mathcal{D})$ is well-defined and that there exists $C>0$ such that,
for any $a,b\geqslant 0$ and any $\mathcal{D}\in\mathfrak{D}_{p}^{a,b}$, one
has

$0\leqslant\Lambda(\mathcal{D})\leqslant C.$

Moreover, for any $1\leqslant p\leqslant\infty$, the subset
$\mathfrak{D}_{p}^{a,b}$ is weakly closed for the norm $L^{p}_{V}$ defined in
(11) (this follows directly from the weak closedness of standard Lebesgue
spaces).

A first useful property to prove the well posedness of the maximization
problem (14) is the following classical result on the concavity of the
functional $\Lambda$, directly obtained from the fact that
$\Lambda(\mathcal{D})$ is an infimum of linear functionals in $\mathcal{D}$
(see for instance [38, Theorem 10.1.1]).

######  Lemma 3.

The functional $\mathcal{D}\mapsto\Lambda(\mathcal{D})$ is concave on the set
of bounded, measurable functions $\mathbb{T}^{d}\to\mathcal{S}_{d}^{+}$.

The following theorem, proved in Section A.1 states the well-posedness of the
optimization problem we consider. Note that we assume that $b>0$, as the
upper-semicontinuity of $\Lambda$, a crucial argument for the result to hold,
is not clear in the case $b=0$.

######  Theorem 4.

Fix $p\in[1,+\infty)$. For any $a\in[0,a_{\mathrm{max}}]$ and $b>0$ such that
$ab\leqslant 1$, there exists a solution
$\mathcal{D}^{\star}\in\mathfrak{D}_{p}^{a,b}$ to (14). In addition, the
optimal diffusion matrix $\mathcal{D}^{\star}$ satisfies the following
properties:

  1. (i)

For any open set $\Omega\subset\mathbb{T}^{d}$, there exists $q\in\Omega$ such
that $\mathcal{D}^{\star}(q)\neq 0$;

  2. (ii)

$\displaystyle\int_{\mathbb{T}^{d}}\left|\mathcal{D}^{\star}(q)\right|_{\mathrm{F}}^{p}\,\mathrm{e}^{-pV(q)}\,dq=1$.

Property (ii) states that, as expected in view of (9), the optimal diffusion
matrix saturates the constraint. Property (i) is useful only for $a=0$. It
guarantees that $\mathcal{D}^{\star}$ cannot be identically 0 on open sets,
but a priori $\mathcal{D}^{\star}(q)$ need not be positive definite, and could
vanish on sets of measure 0. If $\mathcal{D}^{\star}$ cancels at some point,
it may be optimal in the sense of our functional $\Lambda$ but not necessarily
in the sense of the spectral gap; see LABEL:rmk:SG_when_\dim_cancels. We
discuss in Section 3.2 the positive definite character of
$\mathcal{D}^{\star}$. In addition, from a numerical perspective, we prove in
Section 4.1.2 that the optimal diffusion matrix is bounded from below by a
positive constant when the minimization is performed over the intersection of
$\mathfrak{D}_{p}^{a,b}$ and a finite dimensional linear space. This lower
bound is numerically observed to decrease in certain cases when the spatial
discretization is refined (see Figure 3 below).

###  3.2 Characterization of uniformly positive definite maximizers
$\mathcal{D}^{\star}$

The objective of this section is to formally characterize the maximizers
$\mathcal{D}^{\star}$ using the Euler–Lagrange condition satisfied by critical
points of the functional $\Lambda$. When the optimal diffusion is uniformly
positive definite, we first show in Section 3.2.1 that
$\Lambda(\mathcal{D}^{\star})$ is degenerate. Therefore, the mapping
$\mathcal{D}\mapsto\Lambda(\mathcal{D})$ is not necessarily differentiable at
$\mathcal{D}^{\star}$. In order to obtain an Euler–Lagrange equation
characterizing $\mathcal{D}^{\star}$, we circumvent this difficulty in Section
3.2.2 by relying on a smooth approximation of the definition of $\Lambda$. The
theoretical results obtained in this section will be corroborated by numerical
experiments in Section 4.

####  3.2.1 The eigenvalue $\Lambda(\mathcal{D}^{\star})$ is degenerate

Let us consider the solution $\mathcal{D}^{\star}$ to the problem (14) for
$p\in(1,+\infty)$,

$a=0,$

and $b>0$. Note that we exclude the case $p=1$ as the map $M\mapsto\left\lvert
M\right\rvert_{\mathrm{F}}$ is not differentiable at any rank-deficient
$M\in\mathcal{M}_{0,b}$. By Theorem 4, there exists a matrix valued function
$\mathcal{D}^{\star}$ which is a maximizer of $\Lambda$ over the constrained
set $\mathfrak{D}_{p}^{0,b}$. In the remainder of this section, we make two s
on $\mathcal{D}^{\star}$:

$\exists
b_{+}>b,\quad\mathcal{D}^{\star}(q)\,\mathrm{e}^{-V(q)}\leqslant\frac{1}{b_{+}}\,\mathrm{Id}_{d}\textrm{
for a.e. }q\in\mathbb{T}^{d}.$ (16)

and

$\exists c>0,\quad\mathcal{D}^{\star}(q)\geqslant c\,\mathrm{Id}_{d}\textrm{
for a.e. }q\in\mathbb{T}^{d}.$ (17)

From (16) and (17), the constraints associated to the lower bound and upper
bound in $\mathfrak{D}_{p}^{0,b}$ are not active at the critical point
$\mathcal{D}^{\star}$. A standard argument based on the Euler–Lagrange
equation satisfied by the minimizer of the Rayleigh quotient in (7) shows that
there exists a function $u_{\mathcal{D}^{\star}}\in H^{1}_{0}(\mu)$ such that
$\|u_{\mathcal{D}^{\star}}\|_{L^{2}(\mu)}=1$ and
$\Lambda(\mathcal{D}^{\star})=\int_{\mathbb{T}^{d}}\nabla
u_{\mathcal{D}^{\star}}(q)^{\top}\mathcal{D}^{\star}(q)\nabla
u_{\mathcal{D}^{\star}}(q)\,\mu(q)\,dq$.

In our numerical experiments, we always observe that (16) is satisfied up to
choosing $b$ sufficiently small. As explained below, (17) is in particular
natural for potentials which are $1/k$-periodic with
$k\in\mathbb{N}\setminus\\{0\\}$, and $k$ sufficiently large to approach the
homogenized limit, for which the optimal diffusion is indeed uniformly
positive definite; see Section 5.

We show in this section that, when (17) holds, $\Lambda(\mathcal{D}^{\star})$
is a degenerate eigenvalue of the operator
$-\mathcal{L}_{\mathcal{D}^{\star}}$ (under some regularity s in the case
$d=1$). The proof proceeds by contradiction and relies on the fact that, when
$\Lambda(\mathcal{D}^{\star})$ is non-degenerate and (17) holds, one can write
the Euler–Lagrange equation satisfied for $\mathcal{D}^{\star}$. This leads to
a characterization of the optimal diffusion matrix which contradicts (17). In
order to state the result, we need to distinguish for technical reasons the
cases $d=1$ and $d>1$.

######  Lemma 5.

Choose the Frobenius norm for the matrix norm
$\left|\cdot\right|_{\mathrm{F}}$, consider $p\in(1,+\infty)$ and $a=0$.
Assume that the maximizer $\mathcal{D}^{\star}$ solution to (14) satisfies
(16). Assume moreover that
$\mathcal{D}^{\star}\in\mathcal{C}^{0}(\mathbb{T},\mathbb{R}_{+})$ when $d=1$.
If (17) holds, then the optimal spectral gap $\Lambda(\mathcal{D}^{\star})$ is
a degenerate eigenvalue of the diffusion operator
$-\mathcal{L}_{\mathcal{D}^{\star}}$.

###### Proof.

Suppose that $\mathcal{D}^{\star}\geqslant c\mathrm{Id}_{d}$ for $c>0$. We
denote by $u_{\mathcal{D}^{\star}}$ an eigenvector satisfying
$-\mathcal{L}_{\mathcal{D}^{\star}}u_{\mathcal{D}^{\star}}=\Lambda(\mathcal{D}^{\star})u_{\mathcal{D}^{\star}}$.
Observe first that, for any $\delta\mathcal{D}\in
L^{\infty}(\mathbb{T}^{d},\mathcal{S}_{d})$, where $\mathcal{S}_{d}$ is the
set of real, symmetric, matrices of size $d\times d$ (note that
$\delta\mathcal{D}$ does not output a positive matrix a priori), the matrix
valued function $\mathcal{D}^{\star}+t\delta\mathcal{D}$ has values in the
space of symmetric positive definite matrices when $|t|$ is sufficiently
small.

The proof proceeds by contradiction. If the eigenvalue
$\Lambda(\mathcal{D}^{\star})$ is non-degenerate, the eigenvalue
$\Lambda(\mathcal{D}^{\star}+t\delta\mathcal{D})$ remains non-degenerate and
isolated for $|t|$ sufficiently small, and it is possible to choose the
eigenvector $u_{\mathcal{D}^{\star}+t\delta\mathcal{D}}$ associated with this
eigenvalue so that
$\|u_{\mathcal{D}^{\star}+t\delta\mathcal{D}}\|_{L^{2}(\mu)}=1$ and the
mappings $t\mapsto u_{\mathcal{D}^{\star}+t\delta\mathcal{D}}$ and
$t\mapsto\Lambda(\mathcal{D}^{\star}+t\delta\mathcal{D})$ are analytic in an
open neighborhood of $t=0$ (see, e.g., Theorem II.6.1 and the discussion in
Section VII.3.1 of [42]). The optimal diffusion matrix $\mathcal{D}^{\star}$
is therefore characterized by the following Euler–Lagrange equation: there
exists a Lagrange multiplier $\gamma\in\mathbb{R}$ such that

$\forall\delta\mathcal{D}\in
L^{\infty}(\mathbb{T}^{d},\mathcal{S}_{d}),\qquad\frac{d}{dt}{\Lambda}(\mathcal{D}^{\star}+t\delta\mathcal{D})\Big{|}_{t=0}+\gamma\frac{d}{dt}\Phi_{p}(\mathcal{D}^{\star}+t\delta\mathcal{D})\Big{|}_{t=0}=0,$
(18)

where $\Phi_{p}$ encodes the normalization constraint (11):

$\Phi_{p}(\mathcal{D})=\int_{\mathbb{T}^{d}}\left|\mathcal{D}(q)\right|_{\mathrm{F}}^{p}\,\mathrm{e}^{-pV(q)}\,dq.$
(19)

Note that there are no constraints in (18) related to positivity or
boundedness because of (16)-(17).

Relying on the choice of the Frobenius norm for
$\left|\cdot\right|_{\mathrm{F}}$ to compute the differential of $\Phi_{p}$,
(18) can be rewritten as: for any $\delta\mathcal{D}\in
L^{\infty}(\mathbb{T}^{d},\mathcal{S}_{d})$,

$\int_{\mathbb{T}^{d}}\delta\mathcal{D}(q):\left(\nabla
u_{\mathcal{D}^{\star}}\otimes\nabla
u_{\mathcal{D}^{\star}}\right)\mu(q)\,dq=p\gamma\int_{\mathbb{T}^{d}}\left|\mathcal{D}^{\star}(q)\right|_{\mathrm{F}}^{p-2}\mathcal{D}^{\star}(q):\delta\mathcal{D}(q)\,\mathrm{e}^{-pV(q)}\,dq,$

where $:$ and $\otimes$ are respectively the double contraction and outer
product operators: for any $M_{1},M_{2}\in\mathbb{R}^{d\times d}$ and
$\xi,\zeta\in\mathbb{R}^{d}$,

$M_{1}:M_{2}=\mathrm{Tr}\left(M_{1}^{T}M_{2}\right)=\sum_{i,j=1}^{d}\left[M_{1}\right]_{i,j}\left[M_{2}\right]_{i,j},\qquad\left[\xi\otimes\zeta\right]_{i,j}=\xi_{i}\zeta_{j}.$

To compute the derivative of
${\Lambda}(\mathcal{D}^{\star}+t\delta\mathcal{D})$ around $t=0$, as in the
Hellmann–Feynman theorem of quantum mechanics, we have used that
$du_{\mathcal{D}^{\star}+t\delta\mathcal{D}}/dt$ is orthogonal to
$u_{\mathcal{D}^{\star}+t\delta\mathcal{D}}$ in $L^{2}(\mu)$ (which follows
from the normalization condition
$\|u_{\mathcal{D}^{\star}+t\delta\mathcal{D}}\|_{L^{2}(\mu)}=1$), and the
equality

$\nabla^{*}(\mathcal{D}^{\star}\nabla
u_{\mathcal{D}^{\star}})=\Lambda(\mathcal{D}^{\star})u_{\mathcal{D}^{\star}},$
(20)

where we recall that $\nabla^{*}$ is the adjoint of $\nabla$ on $L^{2}(\mu)$.
We therefore deduce from the Euler–Lagrange condition (18) that the optimal
diffusion matrix $\mathcal{D}^{\star}$ satisfies

$\mathcal{D}^{\star}(q)=\gamma_{p}\left|\mathcal{D}^{\star}(q)\right|_{\mathrm{F}}^{2-p}\mathrm{e}^{(p-1)V(q)}\nabla
u_{\mathcal{D}^{\star}}(q)\otimes\nabla u_{\mathcal{D}^{\star}}(q),$ (21)

where the constant $\gamma_{p}>0$ is determined by the constraint
$\Phi_{p}(\mathcal{D}^{\star})=1$.

When $d\geqslant 2$, we immediately obtain a contradiction with the
$\mathcal{D}^{\star}\geqslant c\mathrm{Id}_{d}$ since (21) implies that
$\mathcal{D}^{\star}(q)$ is rank deficient (in fact, rank 1) for a.e.
$q\in\mathbb{T}^{d}$.

When $d=1$, (21) simplifies (after some straightforward manipulations) as

$\mathcal{D}^{\star}(q)=\widetilde{\gamma}_{p}\mathrm{e}^{V(q)}\left|u^{\prime}_{\mathcal{D}^{\star}}(q)\right|^{2/(p-1)},$
(22)

where the constant $\widetilde{\gamma}_{p}>0$ is determined by
$\Phi_{p}(\mathcal{D}^{\star})=1$. We now show that
$u^{\prime}_{\mathcal{D}^{\star}}$ cancels on $\mathbb{T}$, thus contradicting
the that $\mathcal{D}^{\star}(q)\geqslant c>0$ for all $q\in\mathbb{T}$ (the
inequality is indeed valid for all $q\in\mathbb{T}$, and not up to a set of
Lebesgue measure 0 since $\mathcal{D}^{\star}$ is assumed to be continuous).
To do so, we start by noticing that the continuity of $\mathcal{D}^{\star}$
implies that the eigenfunction $u_{\mathcal{D}^{\star}}$ belongs to
$\mathcal{C}^{1}(\mathbb{T},\mathbb{R})$ (see [38, Remark 1.2.11] and [31,
Theorem 9.15]), so that $u_{\mathcal{D}^{\star}}^{\prime}$ is continuous on
$\mathbb{T}$. Since $u_{\mathcal{D}^{\star}}$ is periodic on $\mathbb{T}$, its
derivative $u_{\mathcal{D}^{\star}}^{\prime}$ cancels at least once on
$\mathbb{T}$. In view of (22), the diffusion $\mathcal{D}^{\star}$ then also
vanishes at this point. The contradiction finally allows to obtain the desired
result for $d=1$ as well. ∎

######  Remark 4.

Let us write more precisely the eigenproblem characterizing the second
eigenpair in the situation where the Euler–Lagrange equation (21) holds, which
is for instance the case when the second eigenvalue is nondegenerate. Since
$\left|\cdot\right|_{\mathrm{F}}$ is the Frobenius norm, (21) implies that
$\left|\mathcal{D}^{\star}(q)\right|_{\mathrm{F}}=\gamma_{p}\mathrm{e}^{V(q)}|\nabla
u_{\mathcal{D}^{\star}}(q)|^{2/(p-1)}$ for some parameter $\gamma_{p}>0$
determined by the normalization condition, so that
$\mathcal{D}^{\star}(q)=\gamma_{p}\mathrm{e}^{V(q)}|\nabla
u_{\mathcal{D}^{\star}}(q)|^{-2(p-2)/(p-1)}\nabla
u_{\mathcal{D}^{\star}}(q)\otimes\nabla u_{\mathcal{D}^{\star}}(q)$. In view
of (20) and the fact that, for any vector valued function $\varphi$, it holds
$\nabla^{*}\varphi=-\mu^{-1}\operatorname{div}(\mu\varphi)$, it then follows
that

$-\gamma_{p}\operatorname{div}\left(|\nabla
u_{\mathcal{D}^{\star}}|^{2/(p-1)}\nabla
u_{\mathcal{D}^{\star}}\right)=\Lambda(\mathcal{D}^{\star})\mathrm{e}^{-V}u_{\mathcal{D}^{\star}}.$
(23)

Up to a weight factor on the right-hand side, Equation (23) is an eigenvalue
problem for the $2p/(p-1)$-Laplacian (see, e.g. [52]). This equation is non-
linear. The existing works on eigenvalue problems for nonlinear operators as
the one appearing in (23) suggest that the solutions should be only
$\mathcal{C}^{1,\alpha}$ (see again [52]). Singularities in the derivatives
are indeed numerically observed in the case of vanishing optimal diffusion
coefficients $\mathcal{D}^{\star}$, see _e.g._ Figures 3(b) and 4(a).

####  3.2.2 Euler–Lagrange equation for $\mathcal{D}^{\star}$

Under (17), the first nonzero eigenvalue for the operator associated with the
maximizer $\mathcal{D}^{\star}$ is degenerate, and one cannot rely on the
Euler–Lagrange equation (21). Yet, a characterization of the optimal diffusion
matrix can still be formally recovered using a smooth-min approximation. This
leads to a generalization of the formula (21), which is corroborated by
numerical experiments, see Section 4.2. The positivity of
$\mathcal{D}^{\star}$ is discussed at the end of this section.

##### Smooth-min approximation of the maximization problem.

Our approach relies on the fact that, for a positive $\alpha$, the map
$\mathsf{m}_{\alpha}:\mathbb{R}^{n}\to\mathbb{R}$ defined by

$\mathsf{m}_{\alpha}(x_{1},\dots,x_{n})=\frac{\displaystyle{\sum_{i=1}^{n}}x_{i}\mathrm{e}^{-\alpha
x_{i}}}{\displaystyle{\sum_{i=1}^{n}}\mathrm{e}^{-\alpha x_{i}}}$

converges to the $\min$ function when $\alpha\to+\infty$, _i.e._

$\forall(x_{1},\dots,x_{n})\in\mathbb{R}^{n},\qquad\mathsf{m}_{\alpha}(x_{1},\dots,x_{n})\xrightarrow[\alpha\to+\infty]{}\min\limits_{1\leqslant
i\leqslant n}x_{i}.$ (24)

We adapt this definition to obtain a smooth approximation of the spectral gap
$\Lambda(\mathcal{D})$, from which one can compute the gradient, and exchange
the derivative/limit when $\alpha\to+\infty$ in order to obtain a formal
characterization of the optimal diffusion matrix similar to (21).

We first construct the counterpart of the map $\mathsf{m}_{\alpha}$ in our
framework. In order to do so, we fix $\mathcal{D}\in
L^{p}_{V}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ such that $\mathcal{D}\geqslant
c\mathrm{Id}_{d}$ for some $c>0$ (see (17)). We denote by $g_{\alpha}\colon
x\mapsto\mathrm{e}^{-\alpha x}$ and $h_{\alpha}\colon x\mapsto
x\mathrm{e}^{-\alpha x}$. By functional calculus, the maps $g_{\alpha}$ and
$h_{\alpha}$ can be generalized to functions of the (positive) self-adjoint
operator $-\mathcal{L}_{\mathcal{D}}$. With an abuse of notation, we write in
this case $g_{\alpha}(\mathcal{D})$ and $h_{\alpha}(\mathcal{D})$ instead of
$g_{\alpha}(-\mathcal{L}_{\mathcal{D}})$ and
$h_{\alpha}(-\mathcal{L}_{\mathcal{D}})$. Since $\mathcal{D}$ is bounded and
uniformly bounded from below by a positive constant (in the sense of symmetric
matrices), the positive operator $-\mathcal{L}_{\mathcal{D}}$ has discrete
spectrum, its minimum eigenvalue being simple and equal to 0, see Section 2.1.
We therefore write

$-\mathcal{L}_{\mathcal{D}}=\sum_{i\geqslant 1}\lambda_{i}P_{i},$ (25)

where $\lambda_{1}=0<\lambda_{2}<\ldots<\lambda_{i}<\ldots$ are the distinct
eigenvalues of $-\mathcal{L}_{\mathcal{D}}$ and $P_{i}$ is the orthogonal
projection onto
$\operatorname{Ker}(-\mathcal{L}_{\mathcal{D}}-\lambda_{i}\mathrm{Id})$. The
multiplicity of $\lambda_{i}$ is denoted by
$N_{i}=\mathrm{dim}\,\mathrm{Ran}(P_{i})\geqslant 1$. For $\alpha>0$, we
define

$f_{\alpha}(\mathcal{D})=\frac{\Tr(h_{\alpha}(\mathcal{D}))}{\Tr(g_{\alpha}(\mathcal{D}))-1}=\frac{\Tr\left((-\mathcal{L}_{\mathcal{D}})\mathrm{e}^{\alpha\mathcal{L}_{\mathcal{D}}}\right)}{\Tr(\mathrm{e}^{\alpha\mathcal{L}_{\mathcal{D}}})-1},$
(26)

where the traces are taken on $L^{2}(\mu)$ (in case the traces are taken on
another space than $L^{2}(\mu)$, this will be made precise with a subscript)
and the operators $g_{\alpha}(\mathcal{D}),h_{\alpha}(\mathcal{D})$ are
defined by functional calculus, _i.e._

$\left\\{\begin{aligned} g_{\alpha}(\mathcal{D})&=\sum_{i\geqslant
1}g_{\alpha}(\lambda_{i})P_{i},\\\ h_{\alpha}(\mathcal{D})&=\sum_{i\geqslant
2}h_{\alpha}(\lambda_{i})P_{i},\end{aligned}\right.$

where we used that $h_{\alpha}(\lambda_{1})=0$. The ratio (26) is well-defined
in view of Lemma 23 in Section B.2 and the fact that
$\Tr\left(g_{\alpha}(\mathcal{D})\right)>1$. Formally, since $\lambda_{1}=0$
and $N_{1}=1$,

$f_{\alpha}(\mathcal{D})=\frac{\displaystyle{\sum_{i\geqslant
2}}N_{i}\lambda_{i}\mathrm{e}^{-\alpha\lambda_{i}}}{\displaystyle{\sum_{i\geqslant
2}}N_{i}\mathrm{e}^{-\alpha\lambda_{i}}},$

so that, again formally (compare with (24)),

$f_{\alpha}(\mathcal{D})=\frac{\lambda_{2}\left(1+\dfrac{1}{N_{2}\lambda_{2}}\displaystyle{\sum_{i\geqslant
3}}N_{i}\lambda_{i}\mathrm{e}^{-\alpha(\lambda_{i}-\lambda_{2})}\right)}{1+\dfrac{1}{N_{2}}\displaystyle{\sum_{i\geqslant
3}}N_{i}\mathrm{e}^{-\alpha(\lambda_{i}-\lambda_{2})}}\xrightarrow[\alpha\to+\infty]{}\lambda_{2}=\min\limits_{i\geqslant
2}\lambda_{i}.$ (27)

This limit is justified after the statement of Lemma 23 in Section B.2 using
the dominated convergence theorem. This shows that the quantity
$f_{\alpha}(\mathcal{D})$ is an approximation of $\Lambda(\mathcal{D})$.

This motivates considering the following approximation of the problem (14):
fix $p\in[1,+\infty)$. For any $a\in[0,a_{\mathrm{max}}]$ and $b>0$ such that
$ab\leqslant 1$,

$\textrm{find }\mathcal{D}^{\star,\alpha}\in\mathfrak{D}_{p}^{a,b}\textrm{
such that
}f_{\alpha}(\mathcal{D}^{\star,\alpha})=\sup_{\mathcal{D}\in\mathfrak{D}_{p}^{a,b}}f_{\alpha}(\mathcal{D}).$
(28)

From (27), we expect that $\mathcal{D}^{\star,\infty}$ (which is the limit of
$\mathcal{D}^{\star,\alpha}$ as $\alpha\to+\infty$ whenever it exists) is the
optimal diffusion $\mathcal{D}^{\star}$ solution to (14): this will be
corroborated by numerical experiments in Section 4.2. One can observe that,
for a fixed $\mathcal{D}$, $f_{\alpha}(\mathcal{D})$ is a convex combination
of the eigenvalues $\lambda_{i}$, which implies that
$f_{\alpha}(\mathcal{D})\geqslant\Lambda(\mathcal{D})$. Therefore, it holds
$f_{\alpha}(\mathcal{D}^{\star,\alpha})\geqslant\Lambda(\mathcal{D}^{\star})$.
One can also show that the map $\alpha\mapsto f_{\alpha}(\mathcal{D})$ is
nonincreasing, either by computing its derivative and using the Cauchy–Schwarz
inequality, or by observing that $f_{\alpha}$ is the opposite of the
derivative of the logarithm of the partition function
$\alpha\mapsto\sum_{i\geqslant 2}N_{i}\mathrm{e}^{-\alpha\lambda_{i}}$, and
using the fact that the second derivative of the logarithm of a partition
function is only composed of a variance term and is therefore nonnegative.
These computations can be made rigorous using the theory developed in Appendix
B.

As in the previous section, we assume that there exists a solution to (28) for
$a=0$ and some $b>0$, which satisfies (see (16)-(17))

$\exists
b_{+}>b,\quad\mathcal{D}^{\star,\alpha}(q)\,\mathrm{e}^{-V(q)}\leqslant\frac{1}{b_{+}}\,\mathrm{Id}_{d}\textrm{
for a.e. }q\in\mathbb{T}^{d},$ (29)

and

$\exists c>0,\quad\mathcal{D}^{\star,\alpha}(q)\geqslant
c\,\mathrm{Id}_{d}\textrm{ for a.e. }q\in\mathbb{T}^{d}.$ (30)

In particular, the constraints associated to the lower bound and the upper
bound in $\mathfrak{D}_{p}^{0,b}$ are not active at the critical point
$\mathcal{D}^{\star,\alpha}$. In our numerical experiments, we observe that
(29) is always satisfied for $b$ sufficiently small. Moreover, we use (30)
below for technical reasons (in order to rigorously define $f_{\alpha}$ and
prove that it is differentiable). However, we will observe numerically that
the characterization that we finally obtain on $\mathcal{D}^{\star,\infty}$
(by sending $\alpha\to+\infty$, see Equation (36) below) seems to also hold
even when (30) is not satisfied.

##### Characterization of the optimal diffusion matrix for a finite value of
$\alpha$.

The map $f_{\alpha}$ defined by (26) is smooth over the set of functions
$\mathcal{D}\in L^{p}_{V}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ such that
$\mathcal{D}\geqslant c\mathrm{Id}_{d}$ for some $c>0$, and therefore its
differential can be computed even when some eigenvalues are degenerate. The
main result of this section is the computation of the differentials of the map
$g_{\alpha}$ and $h_{\alpha}$, whose proof is postponed to Section B.1. Note
that the derivatives of the maps $g_{\alpha}$ and $h_{\alpha}$, namely the
maps $g^{\prime}_{\alpha}$ and $h^{\prime}_{\alpha}$ defined by

$\forall x\in\mathbb{R},\forall\alpha>0,\qquad
g_{\alpha}^{\prime}(x)=-\alpha\mathrm{e}^{-\alpha x},\qquad
h_{\alpha}^{\prime}(x)=(1-\alpha x)\mathrm{e}^{-\alpha x},$

can also be generalized to functions of the self-adjoint operator
$-\mathcal{L}_{\mathcal{D}}$, and we then write
$g_{\alpha}^{\prime}(\mathcal{D})=g_{\alpha}^{\prime}(-\mathcal{L}_{\mathcal{D}})$
and
$h_{\alpha}^{\prime}(\mathcal{D})=h_{\alpha}^{\prime}(-\mathcal{L}_{\mathcal{D}})$.

######  Proposition 6.

Fix $\alpha>0$ and $\mathcal{D}\in
L^{p}_{V}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ such that $\mathcal{D}\geqslant
c\mathrm{Id}_{d}$ for some $c>0$. Then there exists $K\in\mathbb{R}_{+}$ such
that, for any $\delta\mathcal{D}\in
L^{\infty}(\mathbb{T}^{d},\mathcal{S}_{d})$ satisfying
$\left\lVert\delta\mathcal{D}\right\rVert_{L^{\infty}(\mathbb{T}^{d})}\leqslant
c/2$, the operators
$g_{\alpha}(\mathcal{D}+\delta\mathcal{D}),g_{\alpha}(\mathcal{D}),g_{\alpha}^{\prime}(\mathcal{D})\mathcal{L}_{\delta\mathcal{D}}$
and the operators
$h_{\alpha}(\mathcal{D}+\delta\mathcal{D}),h_{\alpha}(\mathcal{D}),h_{\alpha}^{\prime}(\mathcal{D})\mathcal{L}_{\delta\mathcal{D}}$
are trace-class on $L^{2}(\mu)$ and

$\left\\{\begin{aligned}
\left\lvert\Tr(g_{\alpha}(\mathcal{D}+\delta\mathcal{D}))-\Tr(g_{\alpha}(\mathcal{D}))-\Tr\left(g_{\alpha}^{\prime}(\mathcal{D})\left(-\mathcal{L}_{\delta\mathcal{D}}\right)\right)\right\rvert&\leqslant
K\left\lVert\delta\mathcal{D}\right\rVert^{2}_{L^{\infty}(\mathbb{T}^{d})},\\\
\left\lvert\Tr(h_{\alpha}(\mathcal{D}+\delta\mathcal{D}))-\Tr(h_{\alpha}(\mathcal{D}))-\Tr\left(h_{\alpha}^{\prime}(\mathcal{D})\left(-\mathcal{L}_{\delta\mathcal{D}}\right)\right)\right\rvert&\leqslant
K\left\lVert\delta\mathcal{D}\right\rVert^{2}_{L^{\infty}(\mathbb{T}^{d})}.\end{aligned}\right.$
(31)

This result immediately implies a differentiability result for the map
$f_{\alpha}$.

######  Corollary 7.

Fix $\alpha>0$ and $\mathcal{D}\in
L^{p}_{V}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ such that $\mathcal{D}\geqslant
c\mathrm{Id}_{d}$ for some $c>0$. Then there exists $K\in\mathbb{R}_{+}$ such
that, for any $\delta\mathcal{D}\in
L^{\infty}(\mathbb{T}^{d},\mathcal{S}_{d})$ satisfying
$\left\lVert\delta\mathcal{D}\right\rVert_{L^{\infty}(\mathbb{T}^{d})}\leqslant
c/2$,

$\left\lvert
f_{\alpha}(\mathcal{D}+\delta\mathcal{D})-f_{\alpha}(\mathcal{D})-\mathcal{F}_{\alpha}(\mathcal{D})(\delta\mathcal{D})\right\rvert\leqslant
K\left\lVert\delta\mathcal{D}\right\rVert^{2}_{L^{\infty}(\mathbb{T}^{d})},$

where

$\mathcal{F}_{\alpha}(\mathcal{D})(\delta\mathcal{D})=\frac{\Tr\left(h^{\prime}_{\alpha}(\mathcal{D})\left(-\mathcal{L}_{\delta\mathcal{D}}\right)\right)\left[\Tr(g_{\alpha}(\mathcal{D}))-1\right]-\Tr(h_{\alpha}(\mathcal{D}))\Tr\left(g^{\prime}_{\alpha}(\mathcal{D})\left(-\mathcal{L}_{\delta\mathcal{D}}\right)\right)}{\left(\Tr(g_{\alpha}(\mathcal{D}))-1\right)^{2}}.$

Straightforward computations lead to

$\mathcal{F}_{\alpha}(\mathcal{D})(\delta\mathcal{D})=\sum\limits_{i\geqslant
1}\left[\frac{\mathcal{G}_{\alpha}(1-\alpha\lambda_{i})+\alpha\mathcal{H}_{\alpha}}{\mathcal{G}_{\alpha}^{2}}\mathrm{e}^{-\alpha\lambda_{i}}\right]\Tr\left(P_{i}\left(-\mathcal{L}_{\delta\mathcal{D}}\right)\right),$

with

$\mathcal{G}_{\alpha}=\sum\limits_{j\geqslant
2}N_{j}\mathrm{e}^{-\alpha\lambda_{j}},\qquad\mathcal{H}_{\alpha}=\sum\limits_{j\geqslant
2}N_{j}\lambda_{j}\mathrm{e}^{-\alpha\lambda_{j}}.$ (32)

We now formally give an explicit formula of the optimal diffusion matrix
$\mathcal{D}^{\star,\alpha}$ solution to (28), following the same computations
as in the proof of Lemma 5. This requires to introduce a Hilbert basis
$(e_{k})_{k\geqslant 1}$ composed of $L^{2}(\mu)$ normalized eigenvectors of
$-\mathcal{L}_{\mathcal{D}}$, namely

$\forall k,\ell\geqslant 1,\qquad-\mathcal{L}_{\mathcal{D}}e_{k}=\text{\char
21\relax}_{k}e_{k},\qquad\left\langle
e_{k},e_{\ell}\right\rangle_{L^{2}(\mu)}=\delta_{k,\ell},$

where $\delta$ is the Kronecker symbol and the eigenvalues $\text{\char
21\relax}_{1}=0\leqslant\text{\char
21\relax}_{2}\leqslant\dots\leqslant\text{\char 21\relax}_{k}\leqslant\dots$
are ordered and counted with their multiplicity:

$\forall i\geqslant 1,\quad\forall
k\in\left\\{1+\sum_{j=1}^{i-1}N_{j},\dots,\sum_{j=1}^{i}N_{j}\right\\},\qquad\text{\char
21\relax}_{k}=\lambda_{i}.$

We rewrite the spectral decomposition (25) of $-\mathcal{L}_{\mathcal{D}}$ in
terms of this given Hilbert basis and with possibly degenerate eigenvalues:
for a given smooth function $\varphi$,

$-\mathcal{L}_{\mathcal{D}}\varphi=\sum_{k\geqslant 1}\text{\char
21\relax}_{k}\left\langle\varphi,e_{k}\right\rangle_{L^{2}(\mu)}e_{k}.$

Exploiting the Euler–Lagrange equation for the objective function
$f_{\alpha}$, namely using that
$\mathcal{F}_{\alpha}(\mathcal{D}^{\star,\alpha})$ is proportional to the
differential of the constraint function (19):

$\forall\delta\mathcal{D}\in
L^{\infty}(\mathbb{T}^{d},\mathcal{S}_{d}),\qquad\mathcal{F}_{\alpha}(\mathcal{D}^{\star,\alpha})(\delta\mathcal{D})\propto\int_{\mathbb{T}^{d}}\left\lvert\mathcal{D}^{\star,\alpha}(q)\right\rvert_{\mathrm{F}}^{p-2}\mathrm{e}^{-pV(q)}\mathcal{D}^{\star,\alpha}(q):\delta\mathcal{D}(q)\,dq.$

Since

$P_{i}=\sum_{k=1+\sum_{l=1}^{i-1}N_{l}}^{\sum_{l=1}^{i}N_{l}}\left\langle
e_{k},\cdot\right\rangle_{L^{2}(\mu)}e_{k}$

and

$\displaystyle\Tr\left(\left\langle
e_{k},\cdot\right\rangle_{L^{2}(\mu)}e_{k}\nabla^{*}\delta\mathcal{D}\nabla\right)$
$\displaystyle=\int_{\mathbb{T}^{d}}e_{k}(q)\left[\nabla^{*}\delta\mathcal{D}(q)\nabla\right]e_{k}(q)\mu(q)\,dq$
$\displaystyle=\int_{\mathbb{T}^{d}}\nabla
e_{k}(q)^{\top}\delta\mathcal{D}(q)\nabla e_{k}(q)\mu(q)\,dq$
$\displaystyle=\int_{\mathbb{T}^{d}}\delta\mathcal{D}(q):\left(\nabla
e_{k}(q)\otimes\nabla e_{k}(q)\right)\mu(q)\,dq,$

one obtains the following characterization of the optimal diffusion matrix for
a fixed value of $\alpha$:

$\mathcal{D}^{\star,\alpha}(q)=\gamma_{\alpha}\left\lvert\mathcal{D}^{\star,\alpha}(q)\right\rvert_{\mathrm{F}}^{2-p}\mathrm{e}^{(p-1)V(q)}\sum_{k\geqslant
2}\left[\frac{\mathcal{G}_{\alpha}(1-\alpha\text{\char
21\relax}_{k,\alpha})+\alpha\mathcal{H}_{\alpha}}{\mathcal{G}_{\alpha}^{2}}\mathrm{e}^{-\alpha\text{\char
21\relax}_{k,\alpha}}\right]\nabla e_{k,\alpha}(q)\otimes\nabla
e_{k,\alpha}(q),$ (33)

with $\gamma_{\alpha}\in\mathbb{R}$ such that
$\Phi_{p}(\mathcal{D}^{\star,\alpha})=1$, and where $(\text{\char
21\relax}_{k,\alpha})_{k\geqslant 1},(e_{k,\alpha})_{k\geqslant 1}$ are the
eigenelements of the operator $-\mathcal{L}_{\mathcal{D}^{\star,\alpha}}$.
Since $\mathcal{D}^{\star,\alpha}$ is a maximizer of $f_{\alpha}$, the
eigenvalues and eigenvectors also depend on $\alpha$, and the definitions of
the quantities $\mathcal{G}_{\alpha}$ and $\mathcal{H}_{\alpha}$ first defined
in (32) and appearing in (33) are changed accordingly:

$\mathcal{G}_{\alpha}=\sum\limits_{j\geqslant 2}\mathrm{e}^{-\alpha\text{\char
21\relax}_{j,\alpha}},\qquad\mathcal{H}_{\alpha}=\sum\limits_{j\geqslant
2}\text{\char 21\relax}_{j,\alpha}\mathrm{e}^{-\alpha\text{\char
21\relax}_{j,\alpha}}.$

Note that we used that $\nabla e_{1,\alpha}=0$ to start the sum at $k=2$ in
(33) since the eigenspace associated to $\text{\char 21\relax}_{1,\alpha}=0$,
which is one-dimensional, is spanned by constant functions. As a consequence,
one obtains in the case $d=1$ (recalling that we assume (30) so that the power
$1/(p-1)$ is well defined)

$\mathcal{D}^{\star,\alpha}(q)=\widetilde{\gamma}_{\alpha}\mathrm{e}^{V(q)}\left(\sum_{k\geqslant
2}\left[\frac{\mathcal{G}_{\alpha}(1-\alpha\text{\char
21\relax}_{k,\alpha})+\alpha\mathcal{H}_{\alpha}}{\mathcal{G}_{\alpha}^{2}}\mathrm{e}^{-\alpha\text{\char
21\relax}_{k,\alpha}}\right]\left\lvert
e_{k,\alpha}^{\prime}(q)\right\rvert^{2}\right)^{1/(p-1)},$ (34)

with $\widetilde{\gamma}_{\alpha}>0$. This formula is corroborated by
numerical examples below, see Figure 2(b).

##### Formal characterization in the limit $\alpha\to+\infty$.

We now formally compute the limit when $\alpha\to+\infty$ in (33). For
notational simplicity, we consider the case when the limiting problem obtained
for $\alpha\to+\infty$ has a leading nonzero eigenvalue with a degeneracy of
at most 2. It is straightforward to adapt the argument to degeneracies of
arbitrary order. We make the following s:

  * •

all the eigenelements converge when $\alpha\to+\infty$: for any $k\geqslant
1$, there exists $\text{\char 21\relax}_{k,\infty}\in\mathbb{R}_{+}$ such that
$\text{\char 21\relax}_{k,\alpha}\xrightarrow[\alpha\to+\infty]{}\text{\char
21\relax}_{k,\infty}$, and there exists $e_{k,\infty}\in H^{1}_{0}(\mu)$ such
that $e_{k,\alpha}\xrightarrow[\alpha\to+\infty]{}e_{k,\infty}$ in
$H^{1}_{0}(\mu)$;

  * •

if $k\geqslant 4$, then $\text{\char 21\relax}_{k,\infty}>\text{\char
21\relax}_{3,\infty}$ (the eigenspace associated to the eigenvalue
$\text{\char 21\relax}_{2,\infty}$ is either 1 or 2 dimensional);

  * •

when $\text{\char 21\relax}_{2,\infty}=\text{\char 21\relax}_{3,\infty}$, we
further assume that the quantity $\alpha(\text{\char
21\relax}_{3,\alpha}-\text{\char 21\relax}_{2,\alpha})$ converges to
$\eta\in\mathbb{R}_{+}$ when $\alpha\to+\infty$.

In the case when $\text{\char 21\relax}_{2,\infty}<\text{\char
21\relax}_{3,\infty}$, the limit when $\alpha\to+\infty$ of the sum on the
right-hand side of (33) can be computed as

$\displaystyle\sum_{k\geqslant
2}\left[\frac{\mathcal{G}_{\alpha}(1-\alpha\text{\char
21\relax}_{k,\alpha})+\alpha\mathcal{H}_{\alpha}}{\mathcal{G}_{\alpha}^{2}}\mathrm{e}^{-\alpha\text{\char
21\relax}_{k,\alpha}}\right]\nabla e_{k,\alpha}\otimes\nabla e_{k,\alpha}$
(35)
$\displaystyle\underset{\alpha\to+\infty}{\sim}\frac{\mathrm{e}^{-\alpha\text{\char
21\relax}_{2,\infty}}(1-\alpha\text{\char
21\relax}_{2,\infty})+\alpha\text{\char
21\relax}_{2,\infty}\mathrm{e}^{-\alpha\text{\char
21\relax}_{2,\infty}}}{\mathrm{e}^{2\alpha\text{\char
21\relax}_{2,\infty}}}\mathrm{e}^{-\alpha\text{\char
21\relax}_{2,\infty}}\nabla e_{2,\infty}\otimes\nabla e_{2,\infty}=\nabla
e_{2,\infty}\otimes\nabla e_{2,\infty}.$

The characterization of the optimal diffusion matrix (33) becomes in the limit
$\alpha\to+\infty$:

$\mathcal{D}^{\star,\infty}(q)=\gamma_{\infty}\left\lvert\mathcal{D}^{\star,\infty}(q)\right\rvert_{\mathrm{F}}^{2-p}\mathrm{e}^{(p-1)V}\nabla
e_{2,\infty}\otimes\nabla e_{2,\infty},$

with $\gamma_{\infty}>0$ such that $\Phi(\mathcal{D}_{p,\infty}^{*})=1$, which
simplifies in the case $d=1$ to

$\mathcal{D}^{\star,\infty}(q)=\widetilde{\gamma}_{\infty}\mathrm{e}^{V(q)}\left\lvert
e_{2,\infty}^{\prime}(q)\right\rvert^{2/(p-1)},$

with $\widetilde{\gamma}_{\infty}>0$. We therefore recover the
characterizations (21) and (22) of Section 3.2.1.

When $\text{\char 21\relax}_{2,\infty}=\text{\char 21\relax}_{3,\infty}$, the
factor in front of $\nabla e_{2,\alpha}\otimes\nabla e_{2,\alpha}$ in (35) is
asymptotically equivalent to

$\displaystyle\frac{\mathrm{e}^{-\alpha\text{\char
21\relax}_{2,\alpha}}}{\left(\mathrm{e}^{-\alpha\text{\char
21\relax}_{2,\alpha}}+\mathrm{e}^{-\alpha\text{\char
21\relax}_{3,\alpha}}\right)^{2}}\left[\left(\mathrm{e}^{-\alpha\text{\char
21\relax}_{2,\alpha}}+\mathrm{e}^{-\alpha\text{\char
21\relax}_{3,\alpha}}\right)\left(1-\alpha\text{\char
21\relax}_{2,\alpha}\right)+\alpha\left(\text{\char
21\relax}_{2,\alpha}\mathrm{e}^{-\alpha\text{\char
21\relax}_{2,\alpha}}+\text{\char
21\relax}_{3,\alpha}\mathrm{e}^{-\alpha\text{\char
21\relax}_{3,\alpha}}\right)\right]$
$\displaystyle\quad=\frac{1}{\left(1+\mathrm{e}^{-\alpha(\text{\char
21\relax}_{3,\alpha}-\text{\char
21\relax}_{2,\alpha})}\right)^{2}}\left[\left(1+\mathrm{e}^{-\alpha(\text{\char
21\relax}_{3,\alpha}-\text{\char
21\relax}_{2,\alpha})}\right)(1-\alpha\text{\char
21\relax}_{2,\alpha})+\alpha\text{\char 21\relax}_{2,\alpha}+\alpha\text{\char
21\relax}_{3,\alpha}\mathrm{e}^{-\alpha(\text{\char
21\relax}_{3,\alpha}-\text{\char 21\relax}_{2,\alpha})}\right]$
$\displaystyle\quad=\frac{1}{\left(1+\mathrm{e}^{-\alpha(\text{\char
21\relax}_{3,\alpha}-\text{\char
21\relax}_{2,\alpha})}\right)^{2}}\left[1+\mathrm{e}^{-\alpha(\text{\char
21\relax}_{3,\alpha}-\text{\char 21\relax}_{2,\alpha})}+\alpha(\text{\char
21\relax}_{3,\alpha}-\text{\char
21\relax}_{2,\alpha})\mathrm{e}^{-\alpha(\text{\char
21\relax}_{3,\alpha}-\text{\char 21\relax}_{2,\alpha})}\right]$
$\displaystyle\quad\xrightarrow[\alpha\to+\infty]{}\frac{1+\mathrm{e}^{-\eta}+\eta\mathrm{e}^{-\eta}}{(1+\mathrm{e}^{-\eta})^{2}},$

while the factor in front of $\nabla e_{3,\alpha}\otimes\nabla e_{3,\alpha}$
is asymptotically equivalent to
$\frac{\mathrm{e}^{-\eta}(1+\mathrm{e}^{-\eta}-\eta)}{(1+\mathrm{e}^{-\eta})^{2}}$.
Therefore, one obtains the following characterization of the optimal diffusion
matrix in the limit $\alpha\to+\infty$:

$\mathcal{D}^{\star,\infty}(q)=\gamma_{\infty}\left\lvert\mathcal{D}^{\star,\infty}(q)\right\rvert_{\mathrm{F}}^{2-p}\mathrm{e}^{(p-1)V}\left[\nabla
e_{2,\infty}\otimes\nabla
e_{2,\infty}+\frac{\mathrm{e}^{-\eta}(1+\mathrm{e}^{-\eta}-\eta)}{1+\mathrm{e}^{-\eta}+\eta\mathrm{e}^{-\eta}}\nabla
e_{3,\infty}\otimes\nabla e_{3,\infty}\right],$ (36)

with $\gamma_{\infty}\in\mathbb{R}$. In the case $d=1$, and when the diffusion
coefficient is positive, the formula (36) simplifies to

$\mathcal{D}^{\star,\infty}(q)=\widetilde{\gamma}_{\infty}\mathrm{e}^{V(q)}\left(\left\lvert
e_{2,\infty}^{\prime}(q)\right\rvert^{2}+\frac{\mathrm{e}^{-\eta}(1+\mathrm{e}^{-\eta}-\eta)}{1+\mathrm{e}^{-\eta}+\eta\mathrm{e}^{-\eta}}\left\lvert
e_{3,\infty}^{\prime}(q)\right\rvert^{2}\right)^{1/(p-1)},$ (37)

with $\widetilde{\gamma}_{\infty}>0$.

Formulas (36)-(37) suggest that, when the first nonzero eigenvalue
$\text{\char 21\relax}_{2,\infty}$ is degenerate, one should take into account
all the eigenvectors associated with it, and reweight their contributions to
obtain the analytical expression of the optimal diffusion matrix. The
correctness of formula (37) is demonstrated on two numerical examples in
Section 4.2, see Figures 1(b) and 2(e) below.

In the case $d=1$, it is clear from (37) that if the value $\eta$ satisfies
$1+\mathrm{e}^{-\eta}-\eta=0$, then the optimal coefficient vanishes where
$e^{\prime}_{2,\infty}$ vanishes. In that case, the (30) is not satisfied.
However, we still believe that the formula (37) holds, from two numerical
experiments in Section 4.2, see Figure 5 below. The particular value of $\eta$
that satisfies $1+\mathrm{e}^{-\eta}-\eta=0$ is

$\eta^{\star}=1+W\left(\frac{1}{\mathrm{e}}\right)\approx 1.27846,$ (38)

where $W$ is the Lambert $W$ function. In fact, a necessary condition
(actually sufficient if $e_{2,\infty}^{\prime}$ and $e_{3,\infty}^{\prime}$ do
not vanish at the same points) for (37) to define a positive definite
diffusion coefficient is $\eta\in[0,\eta^{\star})$. In practice, we only
observed numerical values of $\eta$ in this range.

##  4 Numerical methods for optimizing the spectral gap

To solve the optimization problems (14) and (28) in practice, we introduce two
finite-dimensional parametrizations: one for the diffusion matrix
$\mathcal{D}$, and one for the eigenfunction $u$. We consider for simplicity
the case of isotropic diffusion matrices

$\mathcal{D}(q)=\mathscr{D}(q)\mathrm{Id}_{d}.$ (39)

In Section 4.1, we introduce the numerical methods w used to solve (14) and
(28). We then present some numerical results obtained with our approach in
Section 4.2.

###  4.1 Discretizations of the optimization problems (14) and (28)

We introduce a piecewise constant parametrization for $\mathscr{D}$ in Section
4.1.1; and next a $\mathbb{P}_{1}$ finite element parametrization for the
eigenfunctions $u$ in Section 4.1.2, in the one dimensional case $d=1$ for
simplicity of exposition. The optimization problem then boils down to a
constrained generalized eigenvalue problem, see Section 4.1.3. We describe in
Section 4.1.4 how to solve the latter problem.

####  4.1.1 Parametrization of the diffusion matrix

In practice, we parametrize the scalar valued function $\mathscr{D}$ in (39)
using a finite-dimensional vector subspace of
$L^{p}(\mathbb{T}^{d},\mathbb{R})$. More precisely, we introduce a set of $N$
non-negative functions $\\{\psi_{1},\ldots,\psi_{N}\\}\subset
L^{p}(\mathbb{T}^{d},\mathbb{R})$, and consider functions of the generic form

$\mathscr{D}_{D}(q)=\sum_{n=1}^{N}D_{n}\psi_{n}(q),\qquad
D=(D_{1},\ldots,D_{N})\in\mathbb{R}^{N}.$ (40)

Nonnegativity conditions for $\mathscr{D}_{D}$ may be cumbersome to write when
the functions $\psi_{n}$ have supports which overlap in a nontrivial way, as
is the case for instance for (tensor products of) trigonometric functions.
This motivates us to choose

$\psi_{n}(q)=\mathbf{1}_{K_{n}}(q),$ (41)

the indicator function of a domain $K_{n}\subset\mathbb{T}^{p}$, with the
condition that $(K_{1},\ldots,K_{N})$ forms a partition of $\mathbb{T}^{p}$
(typically, the sets $K_{n}$ are rectangles obtained from a Cartesian mesh).
In order for the function (40) to be in the set (12) for the choice (41), the
components $(D_{n})_{1\leqslant n\leqslant N}$ should belong to the set

$\delta_{p}^{a,b}=\left\\{D\in\left[a\mathrm{e}^{V_{+}(K_{1})},\frac{1}{b}\mathrm{e}^{V_{-}(K_{1})}\right]\times\dots\times\left[a\mathrm{e}^{V_{+}(K_{N})},\frac{1}{b}\mathrm{e}^{V_{-}(K_{N})}\right]\,\middle|\,\sum_{n=1}^{N}\omega_{p,n}D_{n}^{p}\leqslant
1\right\\},$ (42)

where

$V_{-}(E)=\inf_{E}V,\qquad
V_{+}(E)=\sup_{E}V,\qquad\omega_{p,n}=\int_{K_{n}}\mathrm{e}^{-pV}\left\lvert\mathrm{Id}_{d}\right\rvert_{\mathrm{F}}^{p}>0.$

We assume in the sequel that the conditions

$ab\leqslant\min_{1\leqslant n\leqslant
N}\mathrm{e}^{[V_{-}(K_{n})-V_{+}(K_{n})]},\qquad
a\leqslant\left(\sum_{n=1}^{N}\omega_{p,n}\mathrm{e}^{pV_{+}(K_{n})}\right)^{-1/p},$
(43)

hold, so that the set $\delta_{p}^{a,b}$ is not empty. These conditions can be
ensured by choosing $a,b\geqslant 0$ sufficiently small.

The maximization problem (14) can then be approximated by the finite
dimensional maximization problem

$\textrm{find }D^{\star}\in\delta_{p}^{a,b}\textrm{ such that
}\Lambda(\mathscr{D}_{D^{\star}}\mathrm{Id}_{d})=\sup_{D\in\delta_{p}^{a,b}}\Lambda(\mathscr{D}_{D}\mathrm{Id}_{d}).$
(44)

The well-posedness of this optimization problem is guaranteed by the following
result.

######  Proposition 8.

Fix $N\geqslant 1$ and $p\in[1,+\infty)$, and consider a partition
$(K_{1},\dots,K_{N})$ of $\mathbb{T}^{d}$ and $a,b\geqslant 0$ such that (43)
holds. There exists a solution $D^{\star}$ to (44).

Note that the value $b=0$ is allowed for $N$ finite, in contrast to Theorem 4,
since the inequality in (42) ensures that all the components
$(D_{n})_{1\leqslant n\leqslant N}$ are bounded from above by
$(\min_{1\leqslant n\leqslant N}\omega_{p,n})^{-1/p}$.

###### Proof.

The function $D\mapsto\Lambda(\mathscr{D}_{D}\mathrm{Id}_{d})$ is upper-
semicontinuous on $\delta_{p}^{a,b}$ by arguments similar to the ones used in
Section A.1 (and whatever the choice of the norm on $\mathbb{R}^{N}$). In
addition, $\delta_{p}^{a,b}$ is closed and bounded, hence compact, which,
combined with the upper-semicontinuity of the function to maximize, guarantees
the existence of a maximum and thus proves the result. ∎

Likewise, the maximization problem (28) can be approximated by the finite
dimensional maximization problem

$\textrm{find }D^{\star,\alpha}\in\delta_{p}^{a,b}\textrm{ such that
}f_{\alpha}(\mathscr{D}_{D^{\star,\alpha}}\mathrm{Id}_{d})=\sup_{D\in\delta_{p}^{a,b}}f_{\alpha}(\mathscr{D}_{D}\mathrm{Id}_{d}).$

The well-posedness of this optimization problem can be obtained by similar
arguments.

Notice that one can prove that any solution $D^{\star}$ of (44) is positive
when the domains $K_{n}$ have nonempty interiors (see Section A.2 for the
proof).

######  Proposition 9.

Fix $N\geqslant 1$ and $p\in[1,+\infty)$, and consider a partition
$(K_{1},\dots,K_{N})$ of $\mathbb{T}^{d}$ for which the domains
$(K_{n})_{1\leqslant n\leqslant N}$ have nonempty interiors. Assume that
$a,b\geqslant 0$ satisfy (43), and denote by
${D}^{\star}=(D^{\star}_{1},\ldots,{D}^{\star}_{N})\in\delta_{p}^{a,b}$ a
solution to the maximization problem (44). Then,

$\forall n\in\\{1,\ldots,N\\},\qquad{D}^{\star}_{n}>0.$

Note however that this result does not give any information on the behavior of
the lower bound on the components of $D^{\star}$ in the limit of a vanishing
mesh size. The numerical results we present (see for instance Figure 3)
indicate that the lower bound can indeed converge to 0 as $N\to+\infty$.

####  4.1.2 Finite element approximation of the eigenvalue problem

We discuss in this section how to numerically approximate
$\Lambda(\mathcal{D})$ for a given diffusion matrix $\mathcal{D}$. For
simplicity of exposition, we restrict the presentation to the one dimensional
case $d=1$, the extension to higher dimensional situations posing no
difficulties in principle (but numerically viable only for $d\leqslant 3$ in
terms of computational costs).

The minimization in (7) is performed over a linear space of dimension $I$
obtained by a $\mathbb{P}_{1}$ finite element discretization over a regular
mesh with nodes $q_{i}=i/I$ for $i\in\\{0,\ldots,I-1\\}$ with $q_{-1}\equiv
q_{I-1}$, $q_{0}\equiv q_{I}$ and $q_{1}\equiv q_{I+1}$ to comply with
periodic boundary conditions. The continuous basis functions
$(\varphi_{i})_{i=1,\dots,I}$ have support on $[q_{i-1},q_{i+1}]$ and are
piecewise affine:

$\varphi_{i}(q)=\left\\{\begin{aligned} I(q-q_{i-1})&\quad\text{if}\
q\in[q_{i-1},q_{i}],\\\ I(q_{i+1}-q)&\quad\text{if}\
q\in[q_{i},q_{i+1}].\end{aligned}\right.$

We denote by $S_{I}={\rm Span}(\varphi_{1},\ldots,\varphi_{I})$, and by

$S_{I,0}=\left\\{u\in
S_{I}\,\middle|\,\int_{\mathbb{T}}u(q)\mu(q)\,dq=0\right\\}.$

We can then introduce the following mapping, which corresponds to a finite
element approximation of $\Lambda(\mathscr{D}_{D})$:

$\lambda_{N,I}(D)=\inf_{u\in
S_{I,0}\setminus\\{0\\}}\frac{\displaystyle\int_{\mathbb{T}}\mathscr{D}_{D}(q)u^{\prime}(q)^{2}\mu(q)\,dq}{\displaystyle\int_{\mathbb{T}}u(q)^{2}\mu(q)\,dq}.$
(45)

Let us recall the standard approach to solve (45). Upon writing
$u=\sum_{i=1}^{I}U_{i}\varphi_{i}$, the optimization problem (45) can be
rewritten as

$\inf_{U\in\mathbb{R}^{I}\setminus\\{0\\}}\frac{U^{\top}A(D)U}{U^{\top}BU},\qquad
A(D)=\sum_{n=1}^{N}D_{n}A_{n},$

where the matrices $A_{1},\dots,A_{N}$ and $B\in\mathbb{R}^{I\times I}$ have
entries

$[A_{n}]_{i,j}=\int_{K_{n}}\varphi_{j}^{\prime}(q)\varphi_{i}^{\prime}(q)\mathrm{e}^{-V(q)}\,dq,\qquad
B_{i,j}=\int_{\mathbb{T}}\varphi_{j}(q)\varphi_{i}(q)\mathrm{e}^{-V(q)}\,dq.$
(46)

The Euler–Lagrange condition for the minimization problem (45) is therefore
the following finite dimensional eigenvalue problem:

$A(D)U=\sigma BU,\qquad U^{\top}BU=1.$ (47)

Since $A(D)$ and $B$ are real symmetric matrices and $B$ is positive definite,
the eigenvalue problem (47) admits
$0\leqslant\sigma_{1}(D)\leqslant\sigma_{2}(D)\leqslant\ldots\leqslant\sigma_{I}(D)$
real and nonnegative eigenvalues (counted without multiplicities), with
associated $B$-orthonormal eigenvectors
$\\{U_{1}(D),U_{2}(D),\ldots,U_{I}(D)\\}$. In addition, $\sigma_{1}(D)=0$ and
$U_{1}(D)$ is proportional to $(1,\dots,1)$.

####  4.1.3 The discrete optimization problems

The optimization problem (14) is finally approximated by replacing
$\Lambda(\mathscr{D}_{D})$ in (44) by $\lambda_{N,I}(D)$: fix
$p\in[1,+\infty)$. For any $a,b\geqslant 0$ such that (43) holds,

$\text{find }D^{\star}\in\delta^{p}_{a,b}\text{ such that
}\sigma_{2}(D^{\star})=\sup_{D\in\delta_{p}^{a,b}}\sigma_{2}(D)=\lambda_{N,I}^{\star},$
(48)

where $\sigma_{2}(D)$ is the first nonzero (second-smallest) eigenvalue of the
problem (47).

The following proposition guarantees the existence of a minimizer to (48),
which is moreover positive under some compatibility on the meshes used to
discretize the diffusion coefficient and the eigenvalue problem.

######  Proposition 10.

Fix $p\in[1,+\infty)$. Consider $a,b\geqslant 0$ such that (43) holds. Then,
there exists $D^{\star}\in\delta_{p}^{a,b}$ such that

$\sigma_{2}(D^{\star})=\sup_{D\in\delta_{p}^{a,b}}\sigma_{2}(D)=\lambda_{N,I}^{\star}>0.$

Assume moreover that, for all $n\in\\{1,\ldots,N\\}$, there exists
$i\in\\{0,\dots,I-1\\}$ such that $[q_{i-1},q_{i+1}]\subset\overline{K_{n}}$.
Then, any maximizer satisfies $\displaystyle\min_{1\leqslant n\leqslant
N}D^{\star}_{n}>0$.

###### Proof.

As already mentioned in the proof of Proposition 8, the set $\delta_{p}^{a,b}$
is compact. Thus, to prove the existence of a solution, it suffices to show
that the application $D\mapsto\sigma_{2}(D)$ is upper-semicontinuous; in fact,
we prove that it is continuous. Indeed, the application
$M\in\mathcal{S}_{I}(\mathbb{R})\mapsto\sigma_{2}(M)$ mapping any symmetric
matrix to its second-smallest eigenvalue is continuous on
$\mathcal{S}_{I}(\mathbb{R})$ (see, e.g., [11, Corollary III.2.6]). On the
other hand, $D\mapsto A(D)$ is also continuous. By composition,
$D\mapsto\sigma_{2}(D)$ is continuous on $\delta_{p}^{a,b}$, which proves the
existence of a maximizer.

In addition, the Poincaré inequality (15) guarantees with the choice
$D=\gamma\mathbbm{1}_{N}$ (where $\mathbbm{1}_{N}$ is the $N$-dimensional
vector whose components are all equal to 1, and $\gamma>0$ is a normalization
constant chosen such that $D\in\delta_{p}^{a,b}$) that

$\lambda_{N,I}^{\star}\geqslant\sigma_{2}(\gamma\mathbbm{1}_{N})\geqslant\frac{\gamma}{C_{\mu}}>0.$
(49)

To prove that $D^{\star}_{n}>0$ for all $n\in\\{1,\ldots,N\\}$, we proceed by
contradiction. Assume that there exists $n\in\\{1,\ldots,N\\}$ such that
$D^{\star}_{n}\leqslant 0$. Consider $i\in\\{0,\ldots,I-1\\}$ such that
$[q_{i-1},q_{i+1}]\subset\overline{K_{n}}$. Setting
$\phi_{i}=\varphi_{i}-\int_{\mathbb{T}}\varphi_{i}\,\mu$ (so that $\phi_{i}\in
S_{I,0}$), we obtain

$\sigma_{2}(D^{\star})\leqslant\frac{\displaystyle\int_{\mathbb{T}}D^{\star}(q)|\nabla\phi_{i}(q)|^{2}\mu(q)\,dq}{\displaystyle\int_{\mathbb{T}}\phi_{i}(q)^{2}\mu(q)\,dq}=D^{\star}_{n}\frac{\displaystyle\int_{K_{n}}|\nabla\phi_{i}(q)|^{2}\mu(q)\,dq}{\displaystyle\int_{\mathbb{T}}\phi_{i}(q)^{2}\mu(q)\,dq}\leqslant
0,$

which contradicts (49) and thus proves the result. ∎

The optimization problem (28) is approximated by: fix $p\in[1,+\infty)$. For
any $a,b\geqslant 0$ such that (43) holds,

$\text{find }D^{\star,\alpha}\in\delta^{p}_{a,b}\text{ such that
}F_{\alpha}(D^{\star,\alpha})=\sup_{D\in\delta_{p}^{a,b}}F_{\alpha}(D)=\lambda_{N,I}^{\star,\alpha},$
(50)

where

$F_{\alpha}(D)=\frac{H_{\alpha}(D)}{G_{\alpha}(D)},$ (51)

and

$G_{\alpha}(D)=\sum_{i=2}^{I}\mathrm{e}^{-\alpha\sigma_{i}(D)},\qquad
H_{\alpha}(D)=\sum_{i=2}^{I}\sigma_{i}(D)\mathrm{e}^{-\alpha\sigma_{i}(D)}.$
(52)

The well-posedness of the optimization problem (50) (without the positivity of
the optimizer) can be obtained by similar arguments. Using similar (and even
simpler) arguments as the ones presented after (28), one can show that

$\forall\alpha>0,\qquad\lambda_{N,I}^{\star,\alpha}\geqslant\lambda_{N,I}^{\star},$

and the map $\alpha\mapsto F_{\alpha}(D)$ is nonincreasing on
$\mathbb{R}_{+}^{*}$. However, the monotonicity of the map
$\alpha\mapsto\sigma_{2}(D^{\star,\alpha})$ is not clear: it is observed
numerically that this map is nondecreasing, see Figures 2(c) and 4(b).
Finally, let $(\alpha_{n})_{n\geqslant 0}$ be a sequence of
$\mathbb{R}_{+}^{*}$ such that $\alpha_{n}\to+\infty$ when $n\to+\infty$. The
sequence of maximizers $(D^{\star,\alpha_{n}})_{n\geqslant 0}$ is bounded on
$\delta_{p}^{a,b}$ so that, up to a subsequence (not indicated in the
notations), there exists $D^{\star,\infty}\in\delta_{p}^{a,b}$ such that

$D^{\star,\alpha_{n}}\xrightarrow[n\to+\infty]{}D^{\star,\infty}.$

By definition, it holds

$F_{\alpha_{n}}(D^{\star,\alpha_{n}})=\lambda_{N,I}^{\star,\alpha_{n}}\geqslant\lambda_{N,I}^{\star}\geqslant\sigma_{2}(D^{\star,\infty}).$

By continuity of the maps $D\mapsto\sigma_{i}(D)$, it is then easy to show
that $F_{\alpha_{n}}(D^{\star,\alpha_{n}})\to\sigma_{2}(D^{\star,\infty})$, so
that $\lambda_{N,I}^{\star}=\sigma_{2}(D^{\star,\infty})$ and
$D^{\star,\infty}$ solves (48). Moreover,
$\lambda_{N,I}^{\star,\alpha}\to\lambda_{N,I}^{\star}$ as $\alpha\to+\infty$.
The latter convergence is illustrated in Section 4.2.

In conclusion, solving (50) for large values of $\alpha$ should yield a tight
upper bound on $\lambda_{N,I}^{\star}$. In practice, solving (50) for small
values of $\alpha$ already yields a good approximation in the situations we
have considered.

####  4.1.4 Practical aspects of the implementation

We discuss in this section how to numerically solve (14), approximated as
(48). The numerical solution of the smooth-min approximation (28),
approximated as (50), is made precise in Appendix C.1.

To solve (48) in practice, we used either (i) Sequential Least Squares
Quadratic Programming algorithm (SLSQP), which operates through linearization
of the optimality conditions (see for instance [14, Section 15.1] for an
introduction to such methods), relying on the implementation available in
SciPy [69]; or (ii) the IPOPT solver [70] within the Julia programming
language, using the JuMP framework [27]. In order to use these methods, we
need to compute the gradient of the objective function and the constraint,
which we discuss next.

A first task is to approximate the integrals for the matrix elements in (46)
and the weight factors $\omega_{p,n}$ in (42) by some quadrature method. In
our numerical results, we set $N=I$ and $K_{n}=[(n-1)/N,n/N)$ for any
$n\in\left\\{1,\dots,N\right\\}$, so that, for example, $[A_{n}]_{i,j}$ can be
approximated by (with periodic boundary conditions for the indices)

$N\mathrm{e}^{-V\left(\frac{n-1}{N}\right)}(\delta_{i,n}\delta_{i,j}+\delta_{i,n-1}\delta_{i,j}-\delta_{i,n}\delta_{i,j+1}-\delta_{i,n-1}\delta_{i,j-1}),$
(53)

where $\delta$ is the Kronecker delta function, and $B_{i,j}$ can be
approximated by

$\mathrm{e}^{-V\left(\frac{i-1}{N}\right)}\frac{4\delta_{i,j}+\delta_{i,j-1}+\delta_{i,j+1}}{6N}.$
(54)

We refer to Section C.2 for more details in order to obtain (53)-(54).
Likewise, the weights $\omega_{p,n}$ can be approximated by

$\frac{1}{N}\mathrm{e}^{-pV\left(\frac{n-1}{N}\right)}.$ (55)

To run the optimization algorithms, we next need the gradient of the
constraint

$D\mapsto\sum_{n=1}^{N}\omega_{p,n}D_{n}^{p},$

which is readily computed. A more subtle point is to compute the gradient of
the target function $\lambda_{N,I}(D)=\sigma_{2}(D)$ with the notation
introduced in Section 4.1.2. It is readily computed assuming that the second
eigenvalue $\sigma_{2}(D)$ is simple as the following classical result shows.
Recall that the matrix $B$ is symmetric positive definite, so that the matrix
$B^{-1/2}$ is well-defined.

######  Lemma 11.

Consider $D^{0}\in\mathbb{R}^{N}$ such that the second-smallest eigenvalue
$\sigma_{2}(D^{0})$ of the symmetric positive matrix
$B^{-1/2}A(D^{0})B^{-1/2}$ is simple. Then, $D\mapsto\sigma_{2}(D)$ is
differentiable around $D^{0}$, and the components of its gradient are

$\forall n\in\\{1,\dots,N\\},\qquad\frac{\partial\sigma_{2}}{\partial
D_{n}}(D^{0})=U(D^{0})^{\top}A_{n}U(D^{0}),$ (56)

where $U(D^{0})$ is an eigenvector normalized as in (47).

###### Proof.

The proof involves manipulations similar to the ones used to prove Lemma 5.
When $\sigma_{2}(D^{0})$ is simple, $\sigma_{2}$ is smooth in an open
neighborhood of $D^{0}$, and $D\mapsto U(D)$ can also be constructed to be
smooth, so that, by taking the partial derivative with respect to $D_{n}$ of
the first condition in (47),

$\frac{\partial A}{\partial D_{n}}(D)U(D)+A(D)\frac{\partial U}{\partial
D_{n}}(D)=\frac{\partial\sigma_{2}}{\partial
D_{n}}(D)BU(D)+\sigma_{2}(D)B\frac{\partial U}{\partial D_{n}}(D).$

The desired result follows by multiplying this equality by $U(D)^{\top}$ on
the left, using the normalization $U(D)^{\top}BU(D)=1$ and that
$U(D)^{\top}A(D)=\sigma_{2}(D)U(D)^{\top}B$ by transposing the relation (47)
as the matrices $A(D)$ and $B$ are real symmetric. ∎

Using the approximation considered in (53) (recalling that $I=N$ in this
case), the formula approximating the gradient (56) reads: for any
$n\in\left\\{1,\dots,N\right\\}$,

$\displaystyle\frac{\partial\sigma_{2}}{\partial D_{n}}(D)$
$\displaystyle\approx
N\mathrm{e}^{-V\left(\frac{n-1}{N}\right)}\begin{pmatrix}U(D)_{n-1},U(D)_{n}\end{pmatrix}^{\top}\begin{pmatrix}1&-1\\\
-1&1\end{pmatrix}\begin{pmatrix}U(D)_{n-1},U(D)_{n}\end{pmatrix}$
$\displaystyle=N\mathrm{e}^{-V\left(\frac{n-1}{N}\right)}\left[U(D)_{n}-U(D)_{n-1}\right]^{2}.$
(57)

######  Remark 5.

Lemma 11 provides a formula for the gradient of $\sigma_{2}$ around points
where this function is differentiable, namely for $D\in\mathbb{R}^{N}$ such
that $\sigma_{2}(D)$ is simple. However, formula (56) makes sense even if
$\sigma_{2}(D)$ is degenerate, although the resulting quantity is not
necessarily the gradient of $\sigma_{2}$ in this situation. In practice, we
use this formula for the gradient for any $D$, including if degeneracy
happens. The downside is that the convergence may be slower around degeneracy
points.

###  4.2 Numerical results

We apply the optimization procedure described in Section 4.1 to one-
dimensional examples on the torus, with a normalization constraint
corresponding to $p=2$ and a mesh characterized by $N=I=1000$ (recall that $N$
is the dimension of the approximation space for the diffusion coefficient
$\mathcal{D}$, and that $I$ is the number of basis functions for the
$\mathbb{P}_{1}$ finite element method). In particular, we choose
$K_{n}=[(n-1)/N,n/N)$ for $n\in\left\\{1,\dots,N\right\\}$. Note that we used
$N=I$ for simplicity in the implementation of $A_{n}$ (see (53)) in order to
compute the gradient of $\sigma_{2}$ (see Lemma 11) even though the in
Proposition 10 to ensure positivity does not hold in that case. Similar
results were obtained using $I=kN$ with $k\geqslant 2$ which is a sufficient
condition for the positivity result to hold. The integrals for the matrix
elements in (46) are approximated by (53)-(54), and the weights $\omega_{p,n}$
are approximated by (55). More details on the numerical approximations used
are given in Appendix C. We relax the upper bound constraint on the diffusion
coefficient by setting $b=0$.

In this section, we name diffusion coefficients either the finite-dimensional
vectors $D\in\mathbb{R}^{N}$ or the functions defined on $\mathbb{T}$ by (40).
In each figure, the target distribution is $\mu$ defined by (1); the optimal
diffusion coefficient $D^{\star}$ is obtained numerically by solving (48),
while the optimal diffusion coefficient $D^{\star,\alpha}$ is obtained by
solving (50); the constant diffusion coefficient is $D_{\rm
cst}=\gamma\mathbbm{1}_{N}$ with

$\gamma=\left(\sum_{n=1}^{N}\omega_{p,n}\right)^{-1/p},$

so that $\sum_{n=1}^{N}\omega_{p,n}D_{\mathrm{cst},n}^{p}=1$; and the
homogenized diffusion coefficient corresponds to the finite-dimensional
representation of the proxy $\mathcal{D}_{\mathrm{hom}}^{\star}$ for the
optimal diffusion coefficient obtained in a periodic homogenization limit,
whose expression is given by (74) (see Section 5 for the presentation of the
homogenized limit):

$\forall n\in\left\\{1,\dots,N\right\\},\qquad
D^{\star}_{\mathrm{hom},n}=\mathcal{D}^{\star}_{\rm
hom}\left(\frac{n-1}{N}\right).$ (58)

The initial condition for the diffusion coefficient in the optimization
procedure solving either (48) or (50) is set to $D_{\mathrm{hom}}^{\star}$.

We denote by $D^{\star,\infty}$ the numerical approximation of (37). This
vector is defined componentwise as follows: for
$n\in\left\\{1,\dots,N\right\\}$,

$D^{\star,\infty}_{n}=\widetilde{\gamma}_{\infty}\mathrm{e}^{V((n-1)/N)}\left(\left\lvert
U_{2,n}-U_{2,n-1}\right\rvert^{2}+\frac{\mathrm{e}^{-\upeta}(1+\mathrm{e}^{-\upeta}-\upeta)}{1+\mathrm{e}^{-\upeta}+\upeta\mathrm{e}^{-\upeta}}\left\lvert
U_{3,n}-U_{3,n-1}\right\rvert^{2}\right)^{1/(p-1)},$ (59)

where $U_{2}$ (respectively $U_{3}$) is a normalized eigenvector associated
with $\sigma_{2}(D^{\star})$ (respectively $\sigma_{3}(D^{\star})$) and
$\widetilde{\gamma}_{\infty}$ is such that
$\sum_{n=1}^{N}\omega_{p,n}(D^{\star,\infty}_{n})^{p}=1$. Note that periodic
boundary conditions are imposed in (59), _e.g._ $U_{2,0}=U_{2,N}$. The scalar
$\upeta$ is an approximation of $\eta$ (see (37)), which is defined by
$\upeta=\lim\limits_{\alpha\to+\infty}\alpha(\sigma_{3}(D^{\star,\alpha})-\sigma_{2}(D^{\star,\alpha}))$,
where $\sigma_{2}(D^{\star,\alpha})=\lambda^{\star,\alpha}_{N,I}$ and
$\sigma_{3}(D^{\star,\alpha})$ are the leading two nonzero eigenvalues using
the maximizer $D^{\star,\alpha}$ solving (50).

We present the numerical results by first considering cases where $a=0$ in
Section 4.2.1, and then situations where the lower bound $a$ is positive (see
Section 4.2.2). The results for $a=0$ are further distinguished depending on
the value of the approximation of $\eta$ among the cases $\eta=0$,
$\eta\in(0,\eta^{\star})$ and $\eta=\eta^{\star}$. We finally draw conjectures
from our numerical observations and discuss convergence issues in Section
4.2.3.

####  4.2.1 Results for $a=0$

##### A case where $\eta=0$.

We plot in Figure 1(a) various diffusion coefficients as well as the target
distribution for a four-well potential. We first observe that the optimal
diffusion coefficient $D^{\star}$ takes larger values in regions where the
target distribution is small. The analytical proxy $D_{\mathrm{hom}}^{\star}$
has the same general shape as $D^{\star}$. In fact, the two coefficients are
almost superimposed. This observation is formalized in Section 5.3 below.

The numerical values of the spectral gaps obtained for $D^{\star}$,
$D_{\mathrm{hom}}^{\star}$ and $D_{\rm cst}$ are respectively 30.24, 30.19 and
14.70. This simple example thus demonstrates the improvement obtained by
optimizing the diffusion coefficient in terms of the spectral gap of the
operator, and thus the enhanced convergence rate to equilibrium (further
demonstrated by Monte Carlo simulations in Section 6). It also confirms the
relevance of the homogenized proxy $\mathcal{D}_{\mathrm{hom}}^{\star}$
derived in Section 5, which can be computed explicitly and yields a spectral
gap comparable to the one obtained through the optimization procedure.

The optimization procedure leads to an optimal diffusion coefficient which is
uniformly lower bounded on the torus. In particular, we checked that the
minimum of the optimal diffusion coefficient does not decrease as the mesh
gets finer. A different behavior is observed in other numerical examples, see
Figures 3(a) and 4(a) below.

The first three nonzero eigenvalues obtained for the optimal diffusion
coefficient are $\sigma_{2}(D^{\star})=\lambda^{\star}_{N,I}=30.23813406$,
$\sigma_{3}(D^{\star})=30.23813407$, and $\sigma_{4}(D^{\star})=86.06$. This
means that the second eigenvalue can be considered degenerate, which is in
accordance with Lemma 5 (we expect that both eigenvalues are equal in the
limit $I,N\to+\infty$). Figure 1(b) displays the optimal diffusion coefficient
$D^{\star}$ and the numerical approximation $D^{\star,\infty}$ of (37) given
by (59). We have observed that the two eigenvalues
$\sigma_{2}(D^{\star,\alpha})$ and $\sigma_{3}(D^{\star,\alpha})$ are
essentially the same (up to numerical precision) even for very small values of
$\alpha$, the values of
$\alpha(\sigma_{3}(D^{\star,\alpha})-\sigma_{2}(D^{\star,\alpha}))$ when
$\alpha\in[1,7]$ is at most of order $2.5\times 10^{-5}$. We therefore
consider that $\upeta=0$. With this choice, the two coefficients $D^{\star}$
and $D^{\star,\infty}$ are superimposed, as predicted by (37). This confirms
the relevance of the smooth-min approach used in Section 3.2.2.

(a) Various diffusion coefficients and target distribution $\mu$. The solid
curve in red is the optimal diffusion coefficient $D^{\star}$.

(b) Comparison with the formal Euler–Lagrange characterizations (37)
(approximated by (59)) for the case $\eta=0$.

Figure 1: Results of the optimization procedure on the one-dimensional four-
well potential $V(q)=\cos(8\pi q)$.

##### A case where $\eta\in(0,\eta^{\star})$.

We plot in Figure 2 the numerical results for a one-well potential using the
value $\alpha=1$. In Figure 2(a), we plot the maximizer $D^{\star,\alpha}$.
The optimal diffusion coefficient $D^{\star}$ obtained by solving (48) is very
similar, see Figures 2(b) and 2(e)). The numerical values of the spectral gaps
obtained for $D^{\star,\alpha}$, $D_{\mathrm{hom}}^{\star}$ and $D_{\rm cst}$
are respectively 36.75, 32.43 and 30.47. Note that the optimal spectral gap
$\lambda^{\star}_{N,I}$ is 36.88, and that
$\lambda^{\star,\alpha}_{I,N}=36.94$. This shows that the map $f_{\alpha}$ is
a good approximation of the objective function $\Lambda$, and that the
diffusion coefficient obtained when solving (50) yields a spectral gap
comparable to the one obtained when solving (48).

The diffusion coefficient $D^{\star,\alpha}$ is uniformly lower bounded on
$\mathbb{T}$. However, Lemma 5 does not apply since $D^{\star,\alpha}$ is not
optimal in the sense of maximizing the map $\Lambda$. Indeed, we check that
the first three nonzero eigenvalues obtained for $D^{\star,\alpha}$ are 36.75,
37.24 and 138.98. However, the optimal diffusion coefficient $D^{\star}$
obtained by solving (48) yields a degenerate eigenvalue (we obtain
$\sigma_{2}(D^{\star})\approx\sigma_{3}(D^{\star})\approx 36.88$) which is in
accordance with the result of Lemma 5. This means that checking for eigenvalue
degeneracy when solving (50) does not yield information on the eigenvalue
degeneracy of the optimal diffusion coefficient $D^{\star}$. We check in
Figure 2(b) that the characterization of the optimal diffusion coefficient
using a finite value of $\alpha$, _i.e._ formula (34), is valid, using similar
approximations as for (59): more precisely, for
$n\in\left\\{1,\dots,N\right\\}$, the $n$-th component of the discrete
approximation of (34) is given by

$\widetilde{\gamma}_{\alpha}\mathrm{e}^{V((n-1)/N)}\left(\sum_{k=2}^{N}\left[\frac{G_{\alpha}(D^{\star,\alpha})(1-\alpha\sigma_{k}(D^{\star,\alpha}))+\alpha
H_{\alpha}(D^{\star,\alpha})}{G_{\alpha}(D^{\star,\alpha})^{2}}\mathrm{e}^{-\alpha\sigma_{k}(D^{\star,\alpha})}\right]\left\lvert
U_{k,\alpha,n}-U_{k,\alpha,n-1}\right\rvert^{2}\right)^{1/(p-1)},$

where $G_{\alpha},H_{\alpha}$ are defined in (52), $(U_{k,\alpha})_{2\leqslant
k\leqslant N}$ are the normalized eigenvectors associated with the eigenvalues
$(\sigma_{k}(D^{\star,\alpha}))_{2\leqslant k\leqslant N}$ and
$\widetilde{\gamma}_{\alpha}$ is a normalizing constant. Note that this
formula requires all the eigenelements which can be computationally heavy.

We again seek to compute the value of $\upeta$ in order to compare the optimal
diffusion with the formal Euler–Lagrange formula (37). We plot in Figure 2(c)
the two leading nonzero eigenvalues using the maximizer $D^{\star,\alpha}$ for
increasing values of $\alpha$. The behavior is very different from the
previous case: the eigenvalues seem to converge to the same value but are
quite different, even for rather large values of $\alpha$. In Figure 2(d), we
plot the behavior of
$\alpha(\sigma_{3}(D^{\star,\alpha})-\sigma_{2}(D^{\star,\alpha}))$ as
$\alpha$ increases, and fit the data with a function $\alpha\mapsto
K/\alpha+\upeta$, parameterized by $K,\upeta\in\mathbb{R}$, using the method
curve_fit from the scipy.optimize module. This gives a limiting value
$\upeta\approx-0.51$. The comparison between $D^{\star}$ and
$D^{\star,\infty}$ (as in (59)) is given in Figure 2(e): the curves are
superimposed, which again confirms that the formal discussion in Section 3.2.2
is relevant.

(a) Various diffusion coefficients and target distribution $\mu$. The solid
curve in red is the optimal diffusion coefficient $D^{\star,\alpha}$ with
$\alpha=1.0$.

(b) Comparison with the Euler–Lagrange equation (34) using $\alpha=1.0$. The
solid curve in red is the optimal diffusion coefficient $D^{\star}$.

(c) Second and third eigenvalues associated with the maximizer
$D^{\star,\alpha}$ as a function of $\alpha$.

(d) $\alpha(\sigma_{3}(D^{\star,\alpha})-\sigma_{2}(D^{\star,\alpha}))$ as a
function of $\alpha$. The fit is given by $K/\alpha+\upeta$ with
$K\approx-0.027$ and $\upeta\approx 0.51$.

(e) Comparison with the formal Euler–Lagrange characterizations (37)
(approximated by (59)) for the case $\eta\in(0,\eta^{\star})$. The solid curve
in red is the optimal diffusion coefficient $D^{\star}$.

Figure 2:  Results of the optimization procedure on the one-dimensional one-
well potential $V(q)=\cos(2\pi q)$.

##### Cases where $\eta=\eta^{\star}$.

We plot in Figure 3(a) the numerical results for a double-well potential. We
observe that the optimal diffusion coefficient displays two singularities at
two points where it (almost) vanishes, around $q\approx 0.37$ and $q\approx
0.88$. We zoom around the point $q\approx 0.88$ in Figure 3(b) for various
meshes and observe that, as the mesh gets finer ($I=N$ increases), the minimum
value of the optimal diffusion coefficient gets smaller. The analytical proxy
$D_{\mathrm{hom}}^{\star}$ has the same general shape as $D^{\star}$ but
displays no singularity, which may have numerical advantages. The numerical
values of the spectral gaps obtained for $D^{\star}$,
$D_{\mathrm{hom}}^{\star}$ and $D_{\rm cst}$ are respectively 11.23, 10.58 and
0.81. Other numerical properties of the optimal diffusion coefficient were
difficult to obtain for this particular example, as the fact that the optimal
diffusion almost vanishes makes the optimization either unstable or harder to
converge, see Remark 6 below.

We present in Figure 4(a) the numerical results obtained for another double-
well potential. As in the previous example, this diffusion coefficient almost
vanishes and two singularities can be observed around $q\approx 0.25$ and
$q\approx 0.75$. The numerical values of the spectral gap obtained for
$D^{\star}$, $D_{\mathrm{hom}}^{\star}$ and $D_{\rm cst}$ are respectively
22.84, 21.18 and 8.46. The numerical convergence was easier to obtain for this
example. In particular, we were able to compute the value of $\upeta$ by
solving (50) for various values of $\alpha$. We present in Figure 4(b) the
behavior of the two leading nonzero eigenvalues for the values of $\alpha$ for
which the optimization problems converged. We also fit the data similarly to
what was done in Figure 2(d) to extract the limiting value $\upeta\approx
1.276$, which is very close to the critical value $\eta^{\star}$ introduced in
Section 3.2.2, see (38). This corroborates our claim around (38) that when
$\eta=\eta^{\star}$, the optimal diffusion vanishes at some points. We further
back this statement by presenting in Figure 5 the results obtained with (37)
(approximated by (59)) and $\eta=\eta^{\star}$ (the value of $\eta^{\star}$ is
approximated using the lambertw function of the scipy.special module). In
particular, we see in Figure 5(c) that the diffusion given by the
Euler–Lagrange formula with $\eta=\eta^{\star}$ indeed vanishes.

Lastly, note that the singularities observed in Figures 3(b) and 4(a) are
consistent with the discussion in Remark 4. Indeed, when $\eta=\eta^{\star}$,
Equation (21) holds, and we thus expect the optimal diffusion to be non-
smooth.

(a) Various diffusion coefficients and target distribution $\mu$. The solid
curve in red is the optimal diffusion coefficient $D^{\star}$.

(b) Zoom on the optimal diffusion coefficient $D^{\star}$ for increasing
values of $I=N$.

Figure 3:  Results of the optimization procedure on the one-dimensional
double-well potential $V(q)=\sin(4\pi q)(2+\sin(2\pi q))$.

(a) Various diffusion coefficients and target distribution $\mu$. The solid
curve in red is the optimal diffusion coefficient $D^{\star}$.

(b) Second and third eigenvalues associated with the maximizer
$D^{\star,\alpha}$ as a function of $\alpha$. The fit for
$\sigma_{2}(D^{\star,\alpha})$ is the constant $K_{2}\approx 22.842$ and the
fit for $\sigma_{3}(D^{\star,\alpha})$ is given by $\upeta/\alpha+K_{3}$ with
$\upeta\approx 1.276$ and $K_{3}\approx 22.843$.

Figure 4: Results of the optimization procedure on the one-dimensional double-
well potential $V(q)=\cos(4\pi q)$.

(a) $V(q)=\sin(4\pi q)(2+\sin(2\pi q))$, using in (37) the value
$\eta=\eta^{\star}$ given by (38).

(b) $V(q)=\cos(4\pi q)$, using in (37) the value $\eta=\eta^{\star}$ given by
(38).

(c) Zoom around $q\approx 0.75$ for $V(q)=\cos(4\pi q)$, displaying the
coefficient given by (37) with $\eta=\eta^{\star}$ and $\eta=\upeta$.

Figure 5:  Comparison with the formal Euler–Lagrange characterizations (37)
(approximated by (59)) for the cases $\eta=\eta^{\star}$. The solid curves in
red are the optimal diffusion coefficients $D^{\star}$.

####  4.2.2 Results for $a>0$

We now present results of the optimization procedure for the double-well
potential energy function $V(q)=\sin(4\pi q)(2+\sin(2\pi q))$ when imposing
various lowers bounds $a\in\left\\{0,0.2,0.4,0.6,0.8,1\right\\}$. In Figure
6(a), we plot the optimal diffusion coefficients. The case $a=0$ has already
been covered in Figure 3. It is the only one where the optimal diffusion
coefficient almost vanishes on $\mathbb{T}$. When $a=1$, the optimal diffusion
coefficient is fully determined by the normalization constraint, and happens
to coincide with the optimal diffusion coefficient in the homogenized limit
$D_{\rm hom}^{\star}$. We highlight in Figure 6(b) how the lower bounds $a$
affects the shape of the optimal diffusion coefficient: we plot the optimal
diffusion coefficients divided by $D_{\rm hom}^{\star}$, and observe that the
pointwise lower bound $a$ is saturated on larger intervals as $a$ increases to
1. We present in Table 1 the values of the spectral gap for the same values of
the lower bound $a$. Note that the spectral gap for $a=0.2$ is almost the same
as the one for $a=0$. This suggests that imposing a (small) positive lower
bound may be useful in practice, especially in view of the interest of using a
positive diffusion in Markov Chain Monte Carlo algorithms. This point is
further discussed below in Section 6.2 (see Figure 10(c)).

(a) Optimal diffusion coefficients for various lower bounds $a$.

(b) Optimal diffusion coefficients normalized by $D_{\hom}^{\star}$.

Figure 6: Results of the optimization procedure for various lower bounds $a$ and $V(q)=\sin(4\pi q)(2+\sin(2\pi q))$. Lower bound $a$ | 0.0 | 0.2 | 0.4 | 0.6 | 0.8 | 1.0
---|---|---|---|---|---|---
Spectral gap | 11.227 | 11.226 | 11.208 | 11.145 | 10.983 | 10.572
Table 1: Spectral gap associated to optimal diffusion coefficients obtained
with the optimization procedure using different lower bounds $a$ in the same
setting as Figure 6.

####  4.2.3 Numerical conjectures and final remarks

In view of the results above, we conjecture that in all situations, in
dimension 1, for the optimal diffusion coefficient:

  * •

The first three nonzero eigenvalues associated with the operator
$-\mathcal{L}_{\mathcal{D}^{\star}}$ satisfy $\text{\char
21\relax}_{2}=\text{\char 21\relax}_{3}<\text{\char 21\relax}_{4}$.

  * •

The independent and normalized eigenvectors $e_{2}$ and $e_{3}$ (associated to
$\text{\char 21\relax}_{2}=\text{\char 21\relax}_{3}$) have derivatives that
do not vanish at the same points (see (37)), which implies that
$\eta\in[0,\eta^{\star}]$, and there are then two situations:

    * –

either $\eta=\eta^{*}$, which is equivalent to the fact that the optimal
diffusion vanishes at some points;

    * –

or $\eta<\eta^{*}$, which is equivalent to the fact that the optimal diffusion
is uniformly positive.

It would be nice to prove such a conjecture and obtain conditions on $V$
allowing to decide which situation is observed in practice. For example, we
will see in the next section that when $V$ is periodic with a period $1/k$
with an integer $k\in\mathbb{N}\setminus\left\\{0\right\\}$, the optimal
diffusion converges rapidly to the optimal diffusion for the homogenized
limit, which does not vanish.

######  Remark 6.

We have checked that all the numerical results reported in this section are
converged up to numerical precision (with a convergence threshold of
$10^{-15}$). For some potentials, we were not able to get converged results,
in the sense that the convergence threshold was not attained before the
maximum number of iterations, set to 1000 in our examples. In some cases, the
algorithm even stopped before the maximum number of iterations, because it was
unable to find a descent direction compatible with the constraints. In
practice, numerical convergence was harder to obtain when the optimal
diffusion coefficient was almost zero at some points.

##  5 Optimal diffusion matrix in the homogenized limit

A drawback of the optimization method described in Section 4 is the
computational burden associated with the calculation of the optimal diffusion
matrix in practice, which requires solving a high-dimensional eigenvalue
problem. To obtain an approximation of the optimal diffusion, which can be
used for instance as a good initial guess to the optimization procedure, we
propose a simple model for the effective diffusion of the optimized overdamped
Langevin dynamics. This model is obtained by studying the asymptotic behavior
of the optimal diffusion matrix in the homogenized limit of periodic
potentials with infinitely small spatial periods. The limit is explicit either
for linear constraints, or for generic $L^{p}_{V}$ constraints in the one
dimensional case $d=1$, and displays a simple dependence on the target density
$\mu$. The form of this dependence is backed up by prior results in the
literature [61].

In order to present our results, we start by introducing in Section 5.1 the
setting of homogenization theory, in particular the notion of $H$-convergence,
and give some useful and standard results, adapted to our framework of
periodic boundary conditions. Next, in Section 5.2, we study the homogenized
limit of the diffusion matrices and then explicitly optimize the spectral gap
of the homogenized limit for one-dimensional systems or linear constraints.
Finally, in Section 5.3, we study the commutation of the homogenization and
optimization procedures, by showing that the diffusion matrix obtained by
optimizing first and then taking the homogenized limit coincides with the
optimal diffusion matrix for the homogenized problem.

###  5.1 Definitions and compactness results

We present in this section compactness results for sequences of matrices
$(\mathcal{A}^{k})_{k\geqslant 1}\subset
L^{\infty}(\mathbb{T}^{p},\mathcal{M}_{a,b})$. Let us emphasize that we
restrict ourselves to $a,b>0$ in all this section. We start by recalling the
notion of $H$-convergence in the setting of $\mathbb{Z}^{d}$-periodic
functions. Compared to the classical $H$-convergence presented in the
literature (see for instance [2, Chapter 1] and references therein, as well as
[13]), the main difference is that we use periodic boundary conditions rather
than Dirichlet boundary conditions.

######  Definition 1 ($H$-convergence).

Fix $a,b>0$, and consider a sequence $(\mathcal{A}^{k})_{k\geqslant 1}\subset
L^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})$. The sequence
$(\mathcal{A}^{k})_{k\geqslant 1}$ is said to $H$-converge to
$\overline{\mathcal{A}}\in L^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})$
(denoted by $\mathcal{A}^{k}\xrightarrow{H}\overline{\mathcal{A}}$) if, for
any $f\in H^{-1}(\mathbb{T}^{d})$ such that $\langle
f,\mathbf{1}\rangle_{H^{-1}(\mathbb{T}^{d}),H^{1}(\mathbb{T}^{d})}=0$, the
sequence $(u^{k})_{k\geqslant 1}\subset H^{1}(\mathbb{T}^{d})$ of solutions to

$\left\\{\begin{aligned} &-\operatorname{div}\left(\mathcal{A}^{k}\nabla
u^{k}\right)=f&\text{ on }\mathbb{T}^{d},\\\
&\int_{\mathbb{T}^{d}}u^{k}(q)\,dq=0,&\end{aligned}\right.$ (60)

satisfies, in the limit $k\to\infty$,

$\left\\{\begin{aligned} u^{k}&\rightharpoonup u\text{ weakly in
}H^{1}(\mathbb{T}^{d}),\\\ \mathcal{A}^{k}\nabla
u^{k}&\rightharpoonup\overline{\mathcal{A}}\nabla u\text{ weakly in
}L^{2}(\mathbb{T}^{d})^{d},\end{aligned}\right.$

where $u\in H^{1}(\mathbb{T}^{d})$ is the solution of the homogenized problem

$\left\\{\begin{aligned}
&-\operatorname{div}\left(\overline{\mathcal{A}}\nabla u\right)=f&\text{ on
}\mathbb{T}^{d},\\\ &\int_{\mathbb{T}^{d}}u(q)\,dq=0.&\end{aligned}\right.$
(61)

Notice that the functions $u^{k}$ in (60) and $u$ in (61) are indeed well-
defined since $\langle
f,\mathbf{1}\rangle_{H^{-1}(\mathbb{T}^{d}),H^{1}(\mathbb{T}^{d})}=0$ (in
application of the Lax–Milgram theorem, see Lemma 27 in Section D.2). The next
theorem, proved for completeness in Section D.1, states a compactness result
for sequences of periodic functions.

######  Theorem 12 (Murat–Tartar).

Fix $a,b>0$, and consider a sequence $(\mathcal{A}^{k})_{k\geqslant 1}\subset
L^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})$. Up to extraction of a
subsequence, $(\mathcal{A}^{k})_{k\geqslant 1}$ $H$-converges to a matrix
$\overline{\mathcal{A}}\in L^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ (non
constant in general).

The following theorem, adapted from [2, Theorem 1.3.18] (and proved for
completeness in Section D.4.1, see Remark 9), states an $H$-convergence result
for sequences of matrices defined through periodic homogenization. It is a
particular case of Theorem 12 specific to periodic functions.

######  Theorem 13.

Fix $a,b>0$, and let $\mathcal{A}\in
L^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})$. For $k\geqslant 1$, consider
the sequence $(\mathcal{A}_{\\#,k})_{k\geqslant 1}$ defined as

$\forall q\in\mathbb{T}^{d},\qquad\mathcal{A}_{\\#,k}(q)=\mathcal{A}(kq).$

The sequence $(\mathcal{A}_{\\#,k})_{k\geqslant 1}$ $H$-converges to the
constant matrix $\overline{\mathcal{A}}\in\mathcal{M}_{a,b}$ with entries

$\forall 1\leqslant i,j\leqslant
d,\qquad\overline{\mathcal{A}}_{ij}=\int_{\mathbb{T}^{d}}(e_{i}+\nabla
w_{i}(q))^{\top}\mathcal{A}(q)(e_{j}+\nabla w_{j}(q))\,dq,$ (62)

where $(e_{1},\ldots,e_{d})$ is the canonical basis of $\mathbb{R}^{d}$ and
$\\{w_{i}\\}_{1\leqslant i\leqslant d}\subset H^{1}(\mathbb{T}^{d})$ is the
family of unique solutions to the problems

$\left\\{\begin{aligned} &-\operatorname{div}\left(\mathcal{A}(e_{i}+\nabla
w_{i})\right)=0,\\\ &\int_{\mathbb{T}^{d}}w_{i}(q)\,dq=0.\end{aligned}\right.$
(63)

The last theorem, proved in Section D.2 by a straightforward adaptation of [2,
Theorem 1.3.16], states that $H$-convergence implies the convergence of the
sequence of spectral gaps.

######  Theorem 14.

Fix $a,b>0$ and let $(\mathcal{A}^{k})_{k\geqslant 1}\subset
L^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ be a sequence $H$-converging to
a homogenized matrix $\overline{\mathcal{A}}\in
L^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})$. Let $(\rho^{k})_{k\geqslant
1}\subset L^{\infty}(\mathbb{T}^{d})$ be a sequence of positive functions
converging weakly-$*$ in $L^{\infty}(\mathbb{T}^{d})$ to a limiting function
$\rho$, and such that there exists $\rho_{-},\rho_{+}>0$ for which
$0<\rho_{-}\leqslant\rho^{k}(q)\leqslant\rho_{+}<+\infty$ for all $k\geqslant
1$ and for almost all $q\in\mathbb{T}^{d}$. Let $\lambda^{k}$ be the smallest
nonzero eigenvalue and $u^{k}$ an associated normalized eigenvector of the
spectral problem

$\left\\{\begin{aligned} &-\operatorname{div}\left(\mathcal{A}^{k}(q)\nabla
u^{k}(q)\right)=\lambda^{k}\rho^{k}(q)u^{k}(q),\\\
&\int_{\mathbb{T}^{d}}u^{k}(q)^{2}\,dq=1.\end{aligned}\right.$

Then,

$\lim_{k\rightarrow+\infty}\lambda^{k}=\overline{\lambda}$

and, up to a subsequence, $(u^{k})_{k\geqslant 1}$ converges weakly in
$H^{1}(\mathbb{T}^{d})$ to $u$, where $(\overline{\lambda},u)$ is an eigenpair
of the homogenized operator:

$\left\\{\begin{aligned}
&-\operatorname{div}\left(\overline{\mathcal{A}}\nabla
u(q)\right)=\overline{\lambda}\rho(q)u(q),\\\
&\int_{\mathbb{T}^{d}}u(q)^{2}\,dq=1,\end{aligned}\right.$

and $\overline{\lambda}$ is the smallest nonzero eigenvalue of this
homogenized eigenvalue problem.

###  5.2 Optimization of the periodic homogenization limit

We still fix $a,b>0$, and apply the results of Section 5.1 to a particular
class of sequences of matrix valued-functions
$(\mathcal{A}_{\\#,k})_{k\geqslant 1}$, defined for
$\mathcal{D}\in\mathfrak{D}_{p}^{a,b}$ by decreasing the period of the
$\mathbb{Z}^{d}$-periodic matrix

$\mathcal{A}(q)=\mathcal{D}(q)\exp(-V(q))$ (64)

to obtain the $(\mathbb{Z}/k)^{d}$-periodic function

$\mathcal{A}_{\\#,k}(q)=\mathcal{D}_{\\#,k}(q)\exp\left(-V_{\\#,k}(q)\right),\qquad\mathcal{D}_{\\#,k}(q)=\mathcal{D}(kq),\qquad
V_{\\#,k}(q)=V(kq).$ (65)

Note that $\mathcal{A}_{\\#,k}(q)\in\mathcal{M}_{a,b}$ for almost every
$q\in\mathbb{T}^{d}$ and for all $k\geqslant 1$. The periodization procedure
(65) allows us to prove in Section 5.2.1 the $H$-convergence of the sequence
of matrices, and the convergence of the associated eigenvalues and
eigenvectors (see Corollary 15). We then provide in Section 5.2.2 explicit
expressions for diffusion matrices $\mathcal{D}$ which maximize the spectral
gap of the homogenized problem; see in particular Theorem 17 for the one
dimensional case $d=1$.

####  5.2.1 Periodic homogenization of the spectral gap problem

We consider the process $(q^{k}_{t})_{t\geqslant 0}$ defined on
$\mathbb{T}^{d}$ by the dynamics

$dq^{k}_{t}=\left[-\mathcal{D}_{\\#,k}(q^{k}_{t})\nabla
V_{\\#,k}(q^{k}_{t})+\operatorname{div}(\mathcal{D}_{\\#,k})(q^{k}_{t})\right]dt+\sqrt{2}\,\mathcal{D}_{\\#,k}^{1/2}(q^{k}_{t})\,dW_{t}.$

The generator associated to the dynamics acts on test functions
$u:\mathbb{T}^{d}\to\mathbb{R}$ as

$\mathcal{L}_{k}u=\exp\left(V_{\\#,k}\right)\operatorname{div}\left[\mathcal{D}_{\\#,k}\exp\left(-V_{\\#,k}\right)\nabla
u\right].$

Since the eigenfunctions associated with the first eigenvalue 0 are constant
functions, the Rayleigh–Ritz principle implies that the spectral gap of
$\mathcal{L}_{k}$ considered as an operator on
$L^{2}(Z^{-1}\mathrm{e}^{-V_{\\#,k}})$ (recall that $Z$ is defined in (1))
writes (recalling the notation (65)):

$\Lambda_{\\#,k}(\mathcal{D})=\min_{u\in
H^{1}(\mathbb{T}^{d})\setminus\\{0\\}}\left\\{\frac{\displaystyle\int_{\mathbb{T}^{d}}\nabla
u(q)^{\top}\mathcal{A}_{\\#,k}(q)\nabla
u(q)\,dq}{\displaystyle\int_{\mathbb{T}^{d}}u^{2}(q)\,\mathrm{e}^{-V_{\\#,k}(q)}\,dq}\
\middle|\ \int_{\mathbb{T}^{d}}u\,\mathrm{e}^{-V_{\\#,k}}=0\right\\}.$ (66)

Let us emphasize that the test functions to find the minimum on the right-hand
side of the previous equality are $\mathbb{Z}^{d}$-periodic, while the
diffusion and the potential are $(\mathbb{Z}/k)^{d}$-periodic. Note also that
the minimization over the functions $u$ is performed on
$H^{1}(\mathbb{T}^{d})$, which is equivalent in fact to the functional space
$H^{1}(Z^{-1}\mathrm{e}^{-V_{\\#,k}})$ since $\mathbb{T}^{d}$ is compact and
$V$ is smooth. This allows to work with a fixed functional space for the test
functions $u$ whatever the value of $k$.

The spectral gap $\Lambda_{\\#,k}(\mathcal{D})$ is the smallest nonzero
solution $\lambda^{k}$ of the following eigenvalue problem posed on
$\mathbb{T}^{d}$: find $(u^{k},\lambda^{k})\in
H^{1}(\mathbb{T}^{d})\setminus\left\\{0\right\\}\times(0,+\infty)$ such that

$-\operatorname{div}\left[\mathcal{D}_{\\#,k}\exp\left(-V_{\\#,k}\right)\nabla
u^{k}\right]=\lambda^{k}\exp\left(-V_{\\#,k}\right)u^{k},\qquad\int_{\mathbb{T}^{d}}u^{k}\mathrm{e}^{-V_{\\#,k}}=0.$

By applying Theorem 14 to the sequence $(\mathcal{A}_{\\#,k})_{k\geqslant 1}$,
we therefore obtain the following result on the convergence of the sequence of
spectral gaps $\left(\Lambda_{\\#,k}(\mathcal{D})\right)_{k\geqslant 1}$ for a
given diffusion $\mathcal{D}\in\mathfrak{D}_{p}^{a,b}$.

######  Corollary 15.

Fix $p\in[1,+\infty)$, $a,b>0$ satisfying (13) and
$\mathcal{D}\in\mathfrak{D}_{p}^{a,b}$. Then,
$\mathcal{A}_{\\#,k}\xrightarrow{H}\overline{\mathcal{A}}$ with
$\overline{\mathcal{A}}$ defined in (62). Defining
$\overline{\mathcal{D}}=Z^{-1}\overline{\mathcal{A}}$, the sequence
$\left(\Lambda_{\\#,k}(\mathcal{D})\right)_{k\geqslant 1}$ converges as
$k\to\infty$ towards

$\Lambda_{\mathrm{hom}}(\mathcal{D}):=\inf_{u\in
H^{1}_{0}(\mathbb{T}^{d})\setminus\\{0\\}}\frac{\displaystyle\int_{\mathbb{T}^{d}}\nabla
u(q)^{\top}\overline{\mathcal{D}}\nabla
u(q)\,dq}{\displaystyle\int_{\mathbb{T}^{d}}u^{2}(q)\,dq}.$

In the next section, we seek to maximize the spectral gap
$\Lambda_{\mathrm{hom}}(\mathcal{D})$ with respect to the baseline diffusion
matrix $\mathcal{D}$ from which $\overline{\mathcal{D}}$ arises as the
$H$-limit of $\mathcal{D}_{\\#,k}(q)\exp(-V_{\\#,k}(q))/Z$. There is of course
a huge body of literature on optimizing homogenized coefficients, for example
in the context of shape optimization of materials, typically by considering
laminated structures which would correspond here to situations where the
diffusion can only take two values; see for instance [65, 35, 9, 4] and [38,
Chapter 10].

It is useful to provide a more explicit expression of the limiting matrix
$\overline{\mathcal{D}}$ in order to maximize
$\Lambda_{\mathrm{hom}}(\mathcal{D})$. Consider
$\xi=\xi_{1}e_{1}+\dots+\xi_{d}e_{d}\in\mathbb{R}^{d}$ and note that

$\xi^{\top}\overline{\mathcal{D}}\xi=Z^{-1}\int_{\mathbb{T}^{d}}(\xi+\nabla
w_{\xi}(q))^{\top}\mathcal{A}(q)(\xi+\nabla w_{\xi}(q))\,dq,$

where $w_{\xi}=\xi_{1}w_{1}+\dots+\xi_{d}w_{d}$ and the family
$(w_{i})_{1\leqslant i\leqslant d}$ is defined in (63). Recalling (64), it
follows that

$\displaystyle\xi^{\top}\overline{\mathcal{D}}\xi$
$\displaystyle=\int_{\mathbb{T}^{d}}\left(\xi+\nabla
w_{\xi}(q)\right)^{\top}\mathcal{D}(q)\left(\xi+\nabla
w_{\xi}(q)\right)\,\mu(q)\,dq$
$\displaystyle=\int_{\mathbb{T}^{d}}\left(\xi^{\top}\mathcal{D}(q)\xi+2\xi^{\top}\mathcal{D}(q)\nabla
w_{\xi}(q)+\nabla w_{\xi}(q)^{\top}\mathcal{D}(q)\nabla
w_{\xi}(q)\right)\mu(q)\,dq.$ (67)

By the weak formulation of (63), it holds, for any $v\in
H^{1}(\mathbb{T}^{d})$,

$\int_{\mathbb{T}^{d}}\xi^{\top}\mathcal{D}(q)\nabla
v(q)\,\mu(q)\,dq=-\int_{\mathbb{T}^{d}}\nabla
w_{\xi}(q)^{\top}\mathcal{D}(q)\nabla v(q)\,\mu(q)\,dq.$ (68)

Upon replacing $v$ by $w_{\xi}$, we obtain with (67) that

$\xi^{\top}\overline{\mathcal{D}}\xi=\xi^{\top}\left(\int_{\mathbb{T}^{d}}\mathcal{D}(q)\,\mu(q)\,dq\right)\xi-\int_{\mathbb{T}^{d}}\nabla
w_{\xi}(q)^{\top}\mathcal{D}(q)\nabla w_{\xi}(q)\,\mu(q)\,dq.$ (69)

####  5.2.2 Optimization of the homogenized limit

We now focus on the optimization with respect to
$\mathcal{D}\in\mathfrak{D}_{p}^{a,b}$ of the homogenized limit
$\Lambda_{\mathrm{hom}}({\mathcal{D}})$ given in Corollary 15, under some
normalization constraint:

$\Lambda_{\mathrm{hom}}^{\star}=\sup_{\mathcal{D}\in\mathfrak{D}_{p}^{a,b}}\Lambda_{\mathrm{hom}}(\mathcal{D}).$
(70)

This optimization problem is well posed by arguments similar to the ones used
in Section 3.1, and using the fact that the mapping
$\mathcal{D}\mapsto\overline{\mathcal{D}}$ in (69) is continuous on
$\mathfrak{D}_{p}^{a,b}$ (relying on the continuity of the mapping
$L^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})\ni\mathcal{A}\mapsto w_{i}\in
H^{1}(\mathbb{T}^{d})$ where $w_{i}$ is the solution of (63)). In the most
general case, one of the difficulties is to understand how the constraint
$\mathcal{D}\in\mathfrak{D}_{p}^{a,b}$ homogenizes into a constraint on
$\overline{\mathcal{D}}$. We start by analyzing a specific example of
constrained set with linear constraints, which is not the constraint
$\|\mathcal{D}\|_{L^{1}_{V}}\leqslant 1$ with the norm defined in (11). The
interest of the linear constraint we consider is that the optimization problem
becomes trivial. We next derive a similar result for the constraint
$\|\mathcal{D}\|_{L^{p}_{V}}\leqslant 1$, but for the one-dimensional case
$d=1$. Both results hold provided $a,b>0$ are sufficiently small so that the
optimal diffusion matrix automatically satisfies the pointwise lower and upper
bounds involving $a$ and $b$.

##### Linear constraint.

We start by considering the following constraint on $\mathcal{D}$, for a given
constant matrix $M\in\mathcal{S}_{d}^{++}$ (with $\mathcal{S}_{d}^{++}$ the
set of real, symmetric, positive, definite matrices of sizes $d\times d$):

$\int_{\mathbb{T}^{d}}\mathcal{D}(q)\,\mu(q)\,dq=M.$ (71)

######  Lemma 16.

Consider $M\in\mathcal{S}_{d}^{++}$, with $M_{-},M_{+}>0$ such that
$M_{-}\mathrm{Id}_{d}\leqslant M\leqslant M_{+}^{-1}\mathrm{Id}_{d}$. Fix
$a,b>0$ such that

$0<a\leqslant M_{-}Z,\qquad 0<b\leqslant\frac{1}{M_{+}Z}.$ (72)

Then, the diffusion $\mathcal{D}^{\star}_{\mathrm{hom}}(q)=M/\mu(q)$ is a
solution to the optimization problem

$\max\left\\{\Lambda_{\rm{hom}}(\mathcal{D})\,\middle|\,\mathcal{D}\in
L_{V}^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b}),\ \
\int_{\mathbb{T}^{d}}\mathcal{D}(q)\,\mu(q)\,dq=M\right\\}.$

###### Proof.

In view of (69), for any $\mathcal{D}\in
L_{V}^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ satisfying (71), it holds

$\forall\xi\in\mathbb{R}^{d},\qquad\xi^{\top}\overline{\mathcal{D}}\xi\leqslant\xi^{\top}\left(\int_{\mathbb{T}^{d}}\mathcal{D}(q)\,\mu(q)\,dq\right)\xi=\xi^{\top}M\xi.$

Thus, for any $\mathcal{D}\in
L_{V}^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ satisfying the constraint
(71),

$\inf_{u\in
H^{1}_{0}(\mathbb{T}^{d})\setminus\\{0\\}}\frac{\displaystyle\int_{\mathbb{T}^{d}}\nabla
u(q)^{\top}\overline{\mathcal{D}}\nabla
u(q)\,dq}{\displaystyle\int_{\mathbb{T}^{d}}u^{2}(q)\,dq}\leqslant\inf_{u\in
H^{1}_{0}(\mathbb{T}^{d})\setminus\\{0\\}}\frac{\displaystyle\int_{\mathbb{T}^{d}}\nabla
u(q)^{\top}M\nabla u(q)\,dq}{\displaystyle\int_{\mathbb{T}^{d}}u^{2}(q)\,dq}.$

Now, note that $\mathcal{D}(q)=M/\mu(q)$ belongs to
$L_{V}^{\infty}(\mathbb{T}^{d},\mathcal{M}_{a,b})$ in view of the choice (72),
and satisfies $\overline{\mathcal{D}}=M$ in view of (69) since $w_{\xi}=0$ for
any $\xi\in\mathbb{R}^{d}$ by (63). This diffusion matrix therefore maximizes
the limiting spectral gap in view of the above inequality, which concludes the
proof. ∎

#####  $L^{p}_{V}$ constraint in the one-dimensional case.

We now restrict ourselves to the case where $d=1$. In this situation,
$\Lambda_{\mathrm{hom}}(\mathcal{D})$ is proportional to the nonnegative real
number $\overline{\mathcal{D}}$. Thus, maximizing
$\Lambda_{\mathrm{hom}}(\mathcal{D})$ boils down to maximizing