# Optimal non-symmetric Fokker-Planck equation for the convergence to a given
equilibrium

Anton Arnold, Beatrice Signorello  Vienna University of Technology, Institute
of Analysis and Scientific Computing, Wiedner Hauptstr. 8-10, A-1040 Wien,
Austria<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

This paper is concerned with finding Fokker-Planck equations in
$\mathbb{R}^{d}$ with the fastest exponential decay towards a given
equilibrium. For a prescribed, anisotropic Gaussian we determine a non-
symmetric Fokker-Planck equation with linear drift that shows the highest
exponential decay rate for the convergence of its solutions towards
equilibrium. At the same time it has to allow for a decay estimate with a
multiplicative constant arbitrarily close to its infimum.

Such an “optimal” Fokker-Planck equation is constructed explicitly with a
diffusion matrix of rank one, hence being hypocoercive. In an
$L^{2}$–analysis, we find that the maximum decay rate equals the maximum
eigenvalue of the inverse covariance matrix, and that the infimum of the
attainable multiplicative constant is 1, corresponding to the high-rotational
limit in the Fokker-Planck drift. This analysis is complemented with numerical
illustrations in 2D, and it includes a case study for time-dependent
coefficient matrices.

KEYWORDS: Fokker-Planck equation, fastest decay, non-symmetric perturbation,
hypocoercivity
MSC 2010: 35Q84, 35B40, 35Q82, 82C31

##  1\. Introduction

The starting point of this paper is a linear, symmetric (or “reversible”)
Fokker-Planck (FP) equation on $\mathbb{R}^{d}$, $d\geq 2$ with a
corresponding, typically anisotropic Gaussian steady state. It is known from
the literature (see [11, 8], e.g.) that the convergence to equilibrium can be
accelerated by adding to the FP-equation non-symmetric perturbations that do
not alter the equilibrium. It is hence a natural goal to find the “optimal
perturbation” (in a sense to be made precise) such that the corresponding
solutions converge the fastest to the fixed steady state. For FP-equations
with fixed or variable diffusion matrices, this problem was studied,
respectively, in [11] and [8]. A closely related question for the (kinetic) 1D
Goldstein-Taylor system was recently studied in [7]: For a fixed (anti-
symmetric) transport operator, the authors found the best (symmetric)
relaxation operator, yielding the fastest exponential decay to equilibrium.
For the same model, but with constant-in-$x$ relaxation, the propagator norm
was previously computed in [9].

While we shall analyze this problem here in a pure PDE context, the origin of
the question comes from a statistical and probabilistic setting: Let a given
potential $V:\mathbb{R}^{d}\rightarrow\mathbb{R}$ satisfy
$\int_{\mathbb{R}^{d}}e^{-V(x)}dx=1$, and define the probability density
function

(1.1) $f_{\infty,V}(x):=e^{-V(x)}.$

To compute expectations with respect to the associated probability measure
$\mu_{V}$, e.g. via a Markov chain Monte Carlo algorithm (see [6]), one needs
to construct an ergodic Markov process $\left(X_{t}\right)_{t\geq 0}$ with the
unique invariant law $\mu_{V}$, i.e.

(1.2) $\operatorname{law}{(X_{t})}\rightarrow\mu_{V},\qquad\mbox{as
}\>t\rightarrow\infty.$

The efficiency of such algorithms can be measured by the speed of convergence
in (1.2). This motivates to pursue the following objective: find the fastest
among all possible processes that sample from the same equilibrium $\mu_{V}$.
A classical way to sample from $\mu_{V}$ is to consider a standard Brownian
motion with drift $-\nabla V$.

The probability density function $f_{t}$ of the process $X_{t}$ at time $t$
then solves the Fokker-Planck equation

(1.3)
$\partial_{t}f_{t}=\operatorname{div}_{x}(\nabla_{x}f_{t}+\nabla_{x}Vf_{t})=:-Lf_{t},\quad
t>0.$

It is symmetric in the sense that its generator $L$ is symmetric in the
Hilbert space

$\mathcal{H}:=L^{2}(\mathbb{R}^{d},f_{\infty,V}^{-1})=\Big{\\{}f:\mathbb{R}^{d}\rightarrow\mathbb{R},\
\text{measurable
s.t.}\int_{\mathbb{R}^{d}}f(x)^{2}f^{-1}_{\infty,V}(x)dx<\infty\Big{\\}}.$

Under appropriate assumptions on the potential $V$ (e.g. if
$\frac{1}{2}|\nabla V(x)|^{2}-\Delta V(x)\rightarrow+\infty$ for
$|x|\rightarrow\infty$, see [13], A.19) it is possible to show that $f_{t}$
converges to the unique equilibrium $f_{\infty,V}$. Moreover, $L$ is coercive
in $V_{0}^{\perp}$ with
$V_{0}:=\mathrm{span}_{\mathbb{R}}{\\{f_{\infty,V}\\}}\subset\mathcal{H}$,
i.e. $\exists\lambda>0$ such that

(1.4) $\langle
Lf,f\rangle_{\mathcal{H}}\geq\lambda||f||_{\mathcal{H}}^{2},\qquad f\in
V_{0}^{\perp}.$

We shall assume in the sequel that this $\lambda$ is chosen as large as
possible, i.e. as the spectral gap of $L$. As a consequence, if $f_{t}$ is a
solution of (1.3), then

(1.5) $||f_{t}-f_{\infty,V}||_{\mathcal{H}}\leq e^{-\lambda
t}||f_{0}-f_{\infty,V}||_{\mathcal{H}},$

for any normalized initial condition $f_{0}\in\mathcal{H}$ (see Proposition 9
in [13]). So we have a purely exponential convergence estimate.

We shall discuss in the next section that it is often possible to improve the
rate of convergence towards $f_{\infty,V}$ by adding a non-reversible
perturbation in (1.3) while preserving the steady state $f_{\infty,V}$ (as
done in [8] and [11]). As a first step we consider the non-symmetric FP-
equation

(1.6)
$\partial_{t}f_{t}=\operatorname{div}_{x}{(\nabla_{x}f_{t}+(\nabla_{x}V+b)f_{t})}=:-L_{b}f_{t},$

with $b=b(x)$ such that $\operatorname{div}_{x}{(be^{-V})}=0$, to keep the
steady state condition $L_{b}f_{\infty,V}=0$ still valid. In this paper we
will only consider Fokker-Planck equations with linear drift, just as in [11,
8]. This corresponds to quadratic potentials

(1.7) $V(x)=\frac{x^{T}K^{-1}x}{2},\ K\in\mathcal{S}^{>0}$

and linear perturbations of the form $b=Ax,\ A\in\mathcal{M}$.

Notation: Here and in the sequel we denote with $\mathcal{M}$ the set of real
$d\times d$ matrices, $\mathcal{S}^{>0}$ (resp. $\mathcal{S}^{\geq 0}$) the
set of positive definite (resp. positive semi-definite) symmetric matrices,
and with $\mathcal{A}$ the set of anti-symmetric matrices. The spectrum of
$A\in\mathcal{M}$ is denoted by $\sigma(A)$. For a symmetric matrix $A$,
$\lambda_{min}(A)$ and $\lambda_{max}(A)$ denote, respectively, its smallest
and largest eigenvalue.

The following lemma (Lemma 1 in [11]) characterizes explicitly the admissible
perturbations in (1.6).

######  Lemma 1.1.

Let $V(x)$ be given by (1.7) and let $b(x)=Ax$ where $A\in\mathcal{M}$. Then

(1.8) $\operatorname{div}_{x}{(be^{-V})}=0\text{ if and only if
}A=JK^{-1}\text{ with some }\ J\in\mathcal{A}.$

Then the non-symmetric Fokker-Planck equation (1.6) becomes

(1.9)
$\partial_{t}f_{t}=\operatorname{div}_{x}{(\nabla_{x}f_{t}+(I_{d}+J)K^{-1}xf_{t})},$

where $I_{d}$ denotes the identity matrix in $\mathcal{M}$, $J\in\mathcal{A}$
is arbitrary, and $f_{\infty,V}$ is still a steady state.

Note that (1.9) still satisfies (1.5) with the same rate $\lambda$ (see §2.4
in [3]), but $\lambda$ may be smaller than the spectral gap of $L_{b}$.
However, the sharp decay rate can be recovered by hypocoercivity tools [13,
5]: Then one finds constants $\tilde{\lambda}>0$ and $c\geq 1$ (depending on
the fixed potential $V$ and the matrix $J$) such that

(1.10) $||f_{t}-f_{\infty,V}||_{\mathcal{H}}\leq
ce^{-\tilde{\lambda}t}||f_{0}-f_{\infty,V}||_{\mathcal{H}},\quad\forall t\geq
0.$

For the reversible FP-dynamics (with $b=0$) the maximal decay rate
$\tilde{\lambda}$ in the estimate (1.10) is $\lambda$, the biggest coercivity
constant in the inequality (1.4). In this case, the multiplicative constant
$c=1$. The advantage of adding a non-reversible perturbation $b$ is to
possibly obtain a larger decay rate $\tilde{\lambda}>\lambda$, at the price of
allowing for a multiplicative constant $c>1$. In fact, the decay rate may be
improved iff $K$ is not a multiple of $I_{d}$, see §3.2, [11].

The question discussed in [11] is the following: Given the potential (1.7),
which is the optimal non-reversible linear FP-equation of the form (1.9) (and
with time-independent coefficients) such that its solutions converge to
$f_{\infty,V}$ with largest decay rate? For the diffusion matrix fixed as
$I_{d}$, as in (1.9), the authors give a complete answer in [11], Theorem 1.
But if one generalizes the question, allowing to vary both $b$ and the
diffusion matrix, as done in [8], the best decay rate from [11] can be
improved further. Finally, one can extend the question further and analyze if
time-dependent coefficients can enhance the decay of linear FP-equations even
more.

Let us put this paper more into context with the literature on entropy methods
and hypocoercivity: The main goal of [1, 2, 3, 4, 5, 13] is to find explicit
and sometimes even optimal decay rates for a given evolution equation. By
contrast, the novelty in [11, 8] and here is to fix an equilibrium density and
then to seek the evolution equation (within a certain class) that yields the
fastest convergence towards the equilibrium.

This paper is organized as follows: In the next section we formulate this
optimization problem and review the results from [8]. In §3 we present the
main result: As the biggest improvement compared to [8], we shall be able to
obtain multiplicative constants in (1.10) that are arbitrarily close to 1. In
§4.1, §4.2 we will elucidate this result on 2D examples, giving sharp decay
estimates and numerical illustrations. Moreover, we identify the non-symmetric
perturbation of the FP-equation as a highly rotating drift term. Then, in §4.3
we discuss the issue of using time-dependent coefficient matrices to
accelerate the decay behavior, mostly focussing on a numerical case study in
2D. Finally, we conclude in §5.

##  2\. Formulation of the optimization problem and existing results

Let $K\in\mathcal{S}^{>0}$ be given. We define the (typically) anisotropic
Gaussian

(2.1)
$f_{\infty,K}(x):=\frac{\operatorname{det}(K)^{-1/2}}{(2\pi)^{d/2}}\mathrm{exp}{\left(-\frac{x^{T}K^{-1}x}{2}\right)},\qquad
x\in\mathbb{R}^{d},$

and the linear Fokker-Planck equation

(2.2)
$\partial_{t}f_{t}=-L_{C,D}f_{t}:=\mathrm{div}_{x}{(D\nabla_{x}f_{t}+Cxf_{t})},\qquad
x\in\mathbb{R}^{d},\ t\in(0,\infty),$

for arbitrary $x$-independent matrices $D\in\mathcal{S}^{\geq 0}$ and
$C\in\mathcal{M}$. Equation (2.2) is a generalization of the non-reversible
(1.9), possibly with a degenerate (i.e. singular) diffusion matrix $D$.
Moreover we define the set

(2.3) $\mathcal{I}(K):=\\{(C,D)\in\mathcal{M}\times\mathcal{S}^{\geq 0},\
\mathrm{Tr}{(D)}\leq d:L_{C,D}f_{\infty,K}=0\\}.$

The next lemma (Lemma 3.1 in [8]; for $D=I_{d}$ also Lemma 1 in [11]) gives a
characterization of the pairs $(C,D)$ in $\mathcal{I}(K)$.

######  Lemma 2.1.

For $K\in\mathcal{S}^{>0}$ fixed, the following two statements are equivalent:

  * •

$(C,D)\in\mathcal{I}(K)$;

  * •

$D\in\mathcal{S}^{\geq 0},\mathrm{Tr}{(D)}\leq d,$ and $\exists
J\in\mathcal{A}$ such that $C=(D+J)K^{-1}.$

In other words, for $K\in\mathcal{S}^{>0}$ given, we have

(2.4) $\mathcal{I}(K)=\\{((D+J)K^{-1},D)\in\mathcal{M}\times\mathcal{S}^{\geq
0}:\ J\in\mathcal{A},\mathrm{Tr}{(D)}\leq d\\},$

and $\mathcal{I}(K)$ is not empty.

Given a fixed covariance matrix $K\in\mathcal{S}^{>0}$ (and hence the fixed
Gaussian $f_{{\infty},K})$, the set $\mathcal{I}(K)$ represents the matrix
pairs $(C,D)$ such that their associated FP-equation admits $f_{\infty,K}$ as
a normalized steady state. But reversely, for such a FP-equation, the
(normalized) steady state $f_{\infty,K}$ does not have to be unique (e.g.
$C=D=\mathrm{diag}{(1,0)}$ admits (2.1) with any
$K=\mathrm{diag}{(1,\kappa)}$, $\kappa>0$). It is known from the literature
(see for example Theorem 3.1, [5]) that the existence of a unique
$L^{1}-$normalized steady state $f_{\infty,K}$ for (2.2) is equivalent to the
following two conditions on $(C,D)\in\mathcal{M}\times\mathcal{S}^{\geq 0}$:

  1. (1)

$C$ is positive stable (i.e., $C$ has a positive spectral gap
$\rho(C):=\min\\{\operatorname{Re}(\lambda)\>:\>\lambda\in\sigma(C)\\}$);

  2. (2)

hypoellipticity of (2.2) (i.e., there is no non-trivial $C^{T}$-invariant
subspace of $\mathrm{ker}{(D)}$).

For our set-up, hypoellipticity can actually be inferred from
$(C,D)\in\mathcal{I}(K)$; and more precisely:

######  Lemma 2.2.

For some fixed $K\in\mathcal{S}^{>0}$, let $(C,D)\in\mathcal{I}(K)$ and
$\rho(C)>0$. Then the corresponding FP-equation (2.2) is hypoelliptic.

###### Proof.

Normalized steady states of (2.2) are Gaussian with its covariance matrix $Q$
satisfying the continuous Lyapunov equation

(2.5) $2D=CQ+QC^{T}\ .$

Since $D\geq 0$ and $\rho(C)>0$, (2.5) has a unique, symmetric and positive
semi-definite solution $Q$ (see, e.g., Theorem 2.2 in [12]), namely $Q=K$.

By the above mentioned equivalence to the uniqueness of the steady state,
(2.2) is hypoelliptic. ∎

For each fixed steady state $f_{\infty,K}$ we now want to answer the following
questions:

  1. (Q1)

Which FP-evolution(s) converge(s) the fastest, i.e. with largest decay rate
$\lambda_{opt}$ to the steady state in the operator norm of $e^{-L_{C,D}t}$ on
$V_{0}^{\perp}\subset\mathcal{H}:=L^{2}(\mathbb{R}^{d},f_{\infty,K}^{-1})$?

  2. (Q2)

Second, when the best decay rate is fixed, what is the infimum of the
multiplicative constant, $c_{inf}$, in the decay estimate (1.10)?

  3. (Q3)

Third, for a fixed $K\in\mathcal{S}^{>0}$ and the corresponding
$\lambda_{opt}$, and for any $c>c_{inf}$, which pair(s) of matrices
$(C_{opt},D_{opt})\in\mathcal{M}\times\in\mathcal{S}^{\geq 0}$ are such that
$e^{-L_{C_{opt},D_{opt}}t}$ yields the convergence estimate (1.10) with the
constants $(\lambda_{opt},c)$?

  4. (Q4)

Forth, for such an optimal pair of matrices, what bound on $C_{opt}$ can be
found, and how does this bound grow w.r.t. to the space dimension $d$?

  5. (Q5)

Could something be gained by allowing for time-dependent matrices $C(t)$,
$D(t)$?

######  Remark 2.3.

We note that, without the additional constraint $\mathrm{Tr}{(D)}\leq d$ in
the definition of $\mathcal{I}(K)$, the problem of finding an optimal
evolution in the above sense would be ill-posed: Indeed, if $f_{t}$ converges
to $f_{\infty,K}$ as $t\to\infty$, then $f^{\alpha}_{t}:=f_{\alpha t}$, for
any $\alpha>0$ and pertaining to $(\alpha C,\alpha D)$, has the same
equilibrium and converges $\alpha$ times faster to it. For this reason, we
shall only consider diffusion matrices with a prescribed bound for the trace,
as in [8]. In probabilistic language it corresponds to the requirement that
the upper bound on the total amount of randomness simultaneously injected in
the system is prescribed, and this bound is equal to the case $D=I_{d}$.

Next we shall optimize the decay rate within the family of FP-equations (2.2)
satisfying $(C,D)\in\mathcal{I}(K)$. But our choice of matrix $C$ will, in
general, differ from the one constructed in [8]. We base this optimization on
the fact that the sharp exponential decay rate of the FP-equation (2.2) equals
$\rho(C)$ (at least for $C$ diagonalizable, see [5], e.g.). Actually, (2.2)
and its associated drift ODE, i.e. $\frac{d}{dt}x=-\widetilde{C}x$, with
$\widetilde{C}:=K^{-1/2}CK^{1/2}$ (and hence $\rho(C)=\rho(\widetilde{C})$),
have an even closer connection, as proven in Theorem 2.3, [4]:

######  Theorem 2.4.

Let $K\in\mathcal{S}^{>0}$ be given. We consider a FP-equation (2.2) with
$(C,D)\in\mathcal{I}(K)$ and $C$ positive stable. Then, the propagator norms
of (2.2) and of its corresponding drift ODE $\frac{d}{dt}x=-\widetilde{C}x$
are equal, i.e.

(2.6)
$\left\|e^{-L_{C,D}t}\right\|_{\mathcal{B}(V_{0}^{\perp})}=\left\|e^{-\widetilde{C}t}\right\|_{\mathcal{B}(\mathbb{R}^{d})},\qquad\forall
t\geq 0,$

where $||\cdot||_{\mathcal{B}(V_{0}^{\perp})}$ denotes the operator norm on
$\mathcal{H}$ and orthogonality is considered w.r.t. $\mathcal{H}$. Moreover,

$||A||_{\mathcal{B}(\mathbb{R}^{d})}:=\sup_{0\neq
x_{0}\in\mathbb{R}^{d}}\frac{\left|\left|Ax_{0}\right|\right|_{2}}{\left|\left|x_{0}\right|\right|_{2}}$

denotes the spectral matrix norm of any matrix $A\in\mathcal{M}$.

This result motivates to investigate the maximum spectral gap of $C$. Indeed,
the next theorem (see Theorem 2.1 in [8]) identifies the maximum spectral gap
of matrices of the form $C=(D+J)K^{-1}$, and its proof (in [8]) provides an
explicit, algorithmic construction of a corresponding, optimal matrix pair
$(C,D)$.

######  Theorem 2.5.

For $K\in\mathcal{S}^{>0}$ given,

(2.7) $\mathrm{max}{\\{\rho(C):\
(C,D)\in\mathcal{I}(K)\\}}=\mathrm{max}(\sigma(K^{-1}))=\mathrm{min}(\sigma(K))^{-1}.$

Concerning the above questions, the article [8] gives the following (partial)
answers: The authors give a complete and positive answer to question Q1,
obtaining the optimal decay rate $\lambda_{opt}=\max(\sigma(K^{-1}))$. Their
optimal pair $(C_{opt},D_{opt})\in\mathcal{I}(K)$ is very degenerate, the rank
of $D_{opt}$ being one (and this will also be the case for our approach
below). But concerning questions Q2 and Q3, they obtain an estimate for the
multiplicative constant that grows dramatically with the dimension (in fact of
order $d^{40d^{2}}$). This is obtained in [8] when considering a FP-equation
with time-independent coefficients, i.e. the equation form introduced in
(2.2). As a remedy, the authors then considered time-dependent coefficients,
using a symmetric FP-equation with the matrices $(K^{-1},I_{d})$ for small
times and a non-symmetric FP-equation for large times. Discontinuous
coefficients were used there for analytical reasons, to improve decay
estimates. But since their estimates are not sharp, it is not clear if time-
dependent coefficients are really able to enhance the decay property of the
exact FP-propagator norm, i.e. the true function of time, without estimates.
We shall return to this question in §4.3 to elucidate question Q5.
While the main result of [8] is presented for the logarithmic relative
entropy, the same argument works also for the $L^{2}$-norm, as already noted
on page 5, [8]:

######  Theorem 2.6 (Theorem 2.2, [8]).

Let $K\in\mathcal{S}^{>0}$ be given.

  1. (a)

For any ${\color[rgb]{0,0,0}\tilde{c}}>1$ it is possible to construct a matrix
pair $(C_{opt},D_{opt})\in\mathcal{I}(K)$ such that, for all normalized
$f_{0}\in\mathcal{H}$ and for all $t_{0}>0$,

(2.8)
$||f_{t}-f_{\infty,K}||^{2}_{\mathcal{H}}\leq{\color[rgb]{0,0,0}\tilde{c}}\frac{\max(\sigma(K^{-1}))}{2t_{0}}e^{-2\max(\sigma(K^{-1}))(t-t_{0})}||f_{0}-f_{\infty}||^{2}_{\mathcal{H}}\,,\qquad
t\geq t_{0},$

where $f_{t}$ solves the following system of FP-equations

(2.9)
$\begin{cases}\partial_{t}f_{t}=\operatorname{div}_{x}(\nabla_{x}f_{t}+K^{-1}xf_{t}),\qquad\qquad\>\>\>\>0\leq
t\leq t_{0}\,,\\\
\partial_{t}f_{t}=\operatorname{div}_{x}(D_{opt}\nabla_{x}f_{t}+C_{opt}xf_{t}),\qquad
t>t_{0}\,.\end{cases}\,$

  2. (b)

For the choice ${\color[rgb]{0,0,0}\tilde{c}}=2$ in part (a), the matrix
$C_{opt}$ can be estimated as

(2.10) $||C_{opt}||_{\mathcal{F}}\leq
4d^{2}\sqrt{\kappa(K)}\,\lambda_{opt}\,,$

where $||\cdot||_{\mathcal{F}}$ denotes the Frobenius norm
$||A||_{\mathcal{F}}:=\sqrt{\mathrm{Tr}\left(A^{T}A\right)}$, and $\kappa(K)$
is the condition number of $K$.

Optimizing the estimate (2.8) w.r.t. the switching time $t_{0}$, and using the
trivial bound $\left\|e^{-L_{C,D}t}\right\|_{\mathcal{B}(V_{0}^{\perp})}\leq
1$ we obtain the following result:

######  Corollary 2.7.

Under the assumptions of Theorem 2.6, and when choosing
$t_{0}:=\min(\sigma(K))/2$, the following estimate holds for all normalized
$f_{0}\in\mathcal{H}$:

(2.11)
$||f_{t}-f_{\infty,K}||^{2}_{\mathcal{H}}\leq||f_{0}-f_{\infty}||^{2}_{\mathcal{H}}\times\begin{cases}1,\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\>\>\>0\leq
t\leq t_{0}\,,\\\
\min\big{\\{}1,\,{\color[rgb]{0,0,0}\tilde{c}}\,\kappa(K)\,e^{1-2\max(\sigma(K^{-1}))t}\big{\\}},\qquad
t>t_{0}\,.\end{cases}$

where $f_{t}$ solves the FP-system (2.9).

Hence, Theorem 2.2 from [8] can only yield multiplicative constants
$c{\color[rgb]{0,0,0}=\sqrt{\tilde{c}\,\kappa(K)\,e}}>\sqrt{\kappa(K)\,e}$,
using the notation of (1.10).

In the next section we shall improve this result in three directions:
Answering question Q2 we shall prove that $c_{inf}$ is always 1, and
concerning question Q3 we shall construct an optimal matrix pair
$(C_{opt}(c),D_{opt}(c))$ for any given $c>1$. Moreover, we shall not need to
split the FP-evolution in time, in contrast to (2.9). Our key ingredient to
obtain an improved result (compared to [11, 8]) is the equality of the
propagator norms of the FP-equation and of its drift ODE, see Theorem 2.4.
This reduces the quest for an optimal decay estimate to an analogous, and
hence easier ODE problem, without having to invoke a hypocoercive entropy
method as in the proof of Theorem 2.2, [8], or the block-diagonal
decomposition of the FP-propagator as in the proof of Proposition 11, [11].
Finally, concerning question Q4 we shall show that our drift matrix
$C_{opt}(c)$ grows like $\mathcal{O}(d^{3/2})$ (for any fixed $c>1$), compared
to an $\mathcal{O}(d^{2})$–growth in [8].

###  2.1. Time-dependent coefficients

In order to analyze also the decay behavior of the split FP-equation (2.9), we
shall next admit in the FP-equation (2.2) time-dependent coefficient matrices:

(2.12)
$\partial_{t}f_{t}=-L(t)f_{t}:=\mathrm{div}_{x}{(D(t)\nabla_{x}f_{t}+C(t)xf_{t})},\qquad
x\in\mathbb{R}^{d},\ t\in(0,\infty).$

Here we assume that each FP-operator $L(t)$, with $t\geq 0$ fixed, admits
$f_{\infty,K}$ as a steady state, and that the covariance matrix
$K\in\mathcal{S}^{>0}$ is given and independent of $t$. Hence, the coefficient
matrices satisfy $(C(t),D(t))\in\mathcal{I}(K)$ $\forall t\geq 0$, and by
Lemma 2.1:

$C(t)=(D(t)+J(t))\,K^{-1},\quad\mbox{with some }J(t)\in\mathcal{A},\>\>\forall
t\geq 0.$

We shall assume $\forall t\geq 0$ that $\rho(C(t))>0$. Hence, by Lemma 2.2
each FP-operator $L(t)$ is hypocoercive. For (2.12), Theorem 2.4 can be
extended: In the following theorem $S(t_{2},t_{1})$ and $T(t_{2},t_{1})$,
$0\leq t_{1}\leq t_{2}<\infty$ will denote, respectively, the propagator
operators for the PDE (2.12) and the ODE (2.13) that map an initial condition
at time $t_{1}$ to the solution at time $t_{2}$.

######  Theorem 2.8.

Let $K\in\mathcal{S}^{>0}$ be given. Let $(C(t),D(t))\in\mathcal{I}(K)$ be
piecewise smooth functions of $t\geq 0$ (where points of discontinuity do not
accumulate), such that the initial value problem for (2.12) admits a unique
solution in $C([0,\infty);\mathcal{H})$. Then, the propagator norms of (2.12)
and of its corresponding drift ODE,

(2.13) $\frac{d}{dt}x=-\widetilde{C}(t)x,\quad t\in(0,\infty),$

where $\widetilde{C}(t):=K^{-1/2}\,C(t)\,K^{1/2}$, are equal, i.e.:

(2.14)
$\left\|S(t_{2},t_{1})\right\|_{\mathcal{B}(V_{0}^{\perp})}=\left\|T(t_{2},t_{1})\right\|_{\mathcal{B}(\mathbb{R}^{d})},\qquad\forall
0\leq t_{1}\leq t_{2}<\infty.$

Since this result is a straightforward extension of Theorem 2.3 in [4], we
shall give only some hints on the notational differences in the Appendix.

##  3\. Main result

The next theorem is the main result of this work. It states the existence of
pairs $(C_{opt},D_{opt})=(C_{opt}(c),D_{opt}(c))\in\mathcal{I}(K)$ that yield
the maximum decay rate of the propagator norm of $e^{-L_{C_{opt},D_{opt}}t}$,
and in parallel yielding a multiplicative constant $c$ arbitrarily close to
$1$.

######  Theorem 3.1.

Let $K\in\mathcal{S}^{>0}$ be given.

  1. (a)

Then, for any constant $c>1$ there exists a pair $(C_{opt},D_{opt})$
$=$$(C_{opt}(c),D_{opt}(c))\in\mathcal{I}(K)$ such that

(3.1)
$\left|\left|e^{-L_{C_{opt},D_{opt}}t}\right|\right|_{\mathcal{B}(V_{0}^{\perp})}\leq
ce^{-\max(\sigma(K^{-1}))t},\qquad t\geq 0.$

  2. (b)

The matrices from part (a) can be estimated as

(3.2)
$\left|\left|C_{opt}\right|\right|_{\mathcal{F}}\leq\lambda_{opt}\,\Big{[}d+\sqrt{\kappa(K)}\,\frac{2\pi
c^{2}}{\sqrt{3}(c^{2}-1)}\,\sqrt{d}\,(d-1)\Big{]},\quad\left|\left|D_{opt}\right|\right|_{\mathcal{F}}=d.$

In the proof we shall build upon the strategy from $\S 3$ in [11], and only
deviate from their strategy in Step 2 below. Nevertheless we outline the full
proof, to make it readable independently.

###### Proof of Theorem 3.1(a)..

We recall that, given any matrix pair $(C,D)$ in ${\mathcal{I}(K)}$, we can
rewrite the drift matrix $C$ (see Lemma 2.1) as

$C=(D+J)K^{-1}=K^{1/2}(\widetilde{D}+\widetilde{J})K^{-1/2},$

where $\widetilde{D}:=K^{-1/2}DK^{-1/2}$ and
$\widetilde{J}=K^{-1/2}JK^{-1/2}$. Moreover it is easy to check that the map
$M\mapsto K^{-1/2}MK^{-1/2}$ is a bijection that leaves $\mathcal{S}^{\geq 0}$
and $\mathcal{A}$ invariant. We split the proof into three steps.

Step 1 We shall construct an optimal pair
($\widetilde{D}_{opt},\widetilde{J}_{opt})$ and investigate the propagator
norm of the ODE-evolution

(3.3) $\frac{d}{dt}x=-\widetilde{C}_{opt}x,\quad
x_{0}:=x(0)\in\mathbb{R}^{d},\ t\geq 0,$

where $\widetilde{C}_{opt}:=\widetilde{D}_{opt}+\widetilde{J}_{opt}$. More
precisely, we shall provide a decay estimate for
$||e^{-\widetilde{C}_{opt}t}||_{\mathcal{B}(\mathbb{R}^{d})}$ by constructing
an appropriate Lyapunov functional (following §2.1 of [1]).

Following the proof of Theorem 2.1 in [8] we recall that $D$ can enable the
maximum decay rate $\lambda_{opt}:=\max{(\sigma(K^{-1}))}$, only if the range
of $D$ is a subset of $\Omega$, i.e. the eigenspace of $K^{-1}$ corresponding
to $\lambda_{opt}$. Hence we let $v\in\mathbb{R}^{d}$ be a normalized
eigenvector of $K^{-1}$ associated to $\lambda_{opt}$. As in [8] we define the
rank-1 matrix $D_{opt}:=d(v\otimes v)\in\mathcal{S}^{\geq 0}$ with
$\mathrm{Tr}(D_{opt})=d$. It follows that

(3.4) $\widetilde{D}_{opt}=dK^{-1/2}(v\otimes
v)K^{-1/2}=d\lambda_{opt}(v\otimes v)=\lambda_{opt}D_{opt},$

and hence

$\frac{\mathrm{Tr}\left(\widetilde{D}_{opt}\right)}{d}=\lambda_{opt}.$

For $\mathrm{dim}(\Omega)>1$, we remark that the choice of
$\widetilde{D}_{opt}$ made in (3.4) is just one simple option, which enables
the decay rate $\lambda_{opt}$.  For the construction of
$\widetilde{J}_{opt}\in\mathcal{A}$ we use a particular basis of
$\mathbb{R}^{d}$: Let $\\{\psi_{k}\\}_{k=1}^{d}$ be an orthonormal basis of
$\mathbb{R}^{d}$ such that the following condition (Lemma 2, [11]) is
satisfied: for all $k\in\\{1,...,d\\}$,

(3.5)
$\langle\psi_{k},\widetilde{D}_{opt}\psi_{k}\rangle=\frac{\mathrm{Tr}\left(\widetilde{D}_{opt}\right)}{d}=\lambda_{opt}.$

The existence of such basis is guaranteed by Proposition 3 in the same paper.
The essence of the basis $\\{\psi_{k}\\}_{k=1}^{d}$ is to provide an
equidistribution of $\mathrm{Tr}\left(\widetilde{D}_{opt}\right)$ into the
directions $\\{\psi_{k}\\}_{k=1}^{d}$, while $\widetilde{D}_{opt}$ has only
rank 1. This is the starting point to enable a uniform (in $x_{0}$ and $t$)
decay estimate of all trajectories of (3.3), see (3.9) below. We observe that
in [11] the hypotheses of Proposition 3 require $\widetilde{D}_{opt}$ to be
invertible. However, this condition can be weakened to
$\widetilde{D}_{opt}\in\mathcal{S}^{\geq 0}$, as already pointed out in [8]:
$\widetilde{D}_{opt}+\epsilon I_{d}\in\mathcal{S}^{>0}$, and
$\epsilon\rightarrow 0^{+}$ yields the above result.

Next, let $0<\lambda_{1}<\cdots<\lambda_{d}$ be arbitrary numbers in
$\mathbb{R}$ that will be chosen later in a suitable way. We define the matrix
$\widetilde{J}_{opt}:=\Psi\widehat{J}_{opt}\Psi^{-1}\in\mathcal{A}$, with
$\Psi:=\left[\psi_{1},...,\psi_{d}\right]$ and $\widehat{J}_{opt}$ is the
anti-symmetric matrix with elements (as in Lemma 2, [11]):

(3.6)
$\left(\widehat{J}_{opt}\right)_{j,k}:=\frac{\lambda_{j}+\lambda_{k}}{\lambda_{j}-\lambda_{k}}\langle\psi_{j},\widetilde{D}_{opt}\psi_{k}\rangle,\qquad\forall
j\neq k,$

and 0 else.

Now, the strategy consists in finding a suitable symmetric matrix
$P\in\mathcal{S}^{>0}$ that defines a modified norm $||\cdot||_{P}$ in
$\mathbb{R}^{d}$ such that the trajectories of the ODE (3.3) decay with pure
exponential decay rate $\lambda_{opt}$ w.r.t. this norm.

Step 2 Let us proceed with the construction of its inverse matrix
$Q:=P^{-1}\in\mathcal{S}^{>0}$. We define $Q:=\Psi\Lambda\Psi^{-1}$, with the
matrix $\Lambda:=\mathrm{diag}{(\lambda_{1},...,\lambda_{d})}$. We observe
that $Q\in\mathcal{S}^{>0}$ due to the orthonormality of $\Psi$ and the
positivity of $\lambda_{i}$. Moreover, by definition, $Q$ has the eigenvectors
$\psi_{i}$ and eigenvalues $\lambda_{i}$. By using Lemma 2, [11] (or a
straightforward computation using (3.5)) the following Lyapunov equation holds
for $Q$:

(3.7)
$\widetilde{J}_{opt}Q-Q\widetilde{J}_{opt}=-Q\widetilde{D}_{opt}-\widetilde{D}_{opt}Q+2\lambda_{opt}Q.$

Let us define the modified norm $||x||^{2}_{P}:=\langle x,Px\rangle$ on
$\mathbb{R}^{d}$, where $P:=Q^{-1}\in\mathcal{S}^{>0}$. Differentiating this
norm along a trajectory of the ODE (3.3) we obtain with (3.7), multiplied on
either side by $P=Q^{-1}$:

(3.8)
$\frac{d}{dt}||x(t)||^{2}_{P}=-\Big{\langle}x(t),\left[P\left(\widetilde{D}_{opt}+\widetilde{J}_{opt}\right)+\left(\widetilde{D}_{opt}-\widetilde{J}_{opt}\right)P\right]x(t)\Big{\rangle}=-2\lambda_{opt}||x(t)||^{2}_{P}.$

Hence the modified norm decays with rate $\lambda_{opt}$, i.e.

(3.9) $||x(t)||_{P}^{2}=e^{-2\lambda_{opt}t}||x(0)||_{P}^{2},\qquad t\geq 0.$

Transforming to the Euclidean vector norm, we obtain for the propagator

(3.10)
$\left|\left|e^{-\widetilde{C}_{opt}t}\right|\right|_{\mathcal{B}(\mathbb{R}^{d})}\leq\sqrt{\kappa(P)}\,e^{-\lambda_{opt}t},\qquad
t\geq 0,$

where $\kappa(P)$ denotes the condition number of the matrix $P$.

Step 3 The multiplicative constant appearing in (3.10) can be adjusted by
choosing the eigenvalues of $P$ in the following way: Given any $c>1$, and due
to the fact that $\kappa(P)=\kappa(Q)=\frac{\lambda_{d}}{\lambda_{1}}$, it is
sufficient to choose $\lambda_{d}$ and $\lambda_{1}$ such that their quotient
is (less or) equal to $c^{2}$. The remaining parameters
$\lambda_{2}<...<\lambda_{d-1}\in(\lambda_{1},\lambda_{d})$ could be freely
chosen at this point, but assigning them a precise value will be crucial in
the proof of part (b).

To summarize, we have proved so far that, for any prescribed $c>1$, there
exists a pair of matrices $\widetilde{J}_{opt}\in\mathcal{A}$ and
$\widetilde{D}_{opt}\in\mathcal{S}^{\geq 0}$ such that

(3.11)
$\left|\left|e^{-\widetilde{C}_{opt}t}\right|\right|_{\mathcal{B}(\mathbb{R}^{d})}\leq
ce^{-\lambda_{opt}t},\qquad t\geq 0.$

We conclude the proof by combining Theorem 2.4 applied to the operator
$e^{-L_{C_{opt},D_{opt}}t}$, and the above inequality (3.11). ∎

For $K=\alpha I_{d}$, we remark that a trivial modification of the above proof
admits the choice $(C_{opt},D_{opt})=(K^{-1},I_{d})$, $J=0$, $P=I_{d}$. In
this case the reversible dynamics is already optimal with
$\lambda_{opt}=\alpha^{-1}$ and $c=1$ in (3.1). Moreover,
$\left|\left|C\right|\right|_{\mathcal{F}}=\lambda_{opt}\sqrt{d}$.

###### Proof of Theorem 3.1(b)..

First we compute the Frobenius norm of $D_{opt}:=d(v\otimes v)$, with
$v\in\mathbb{R}^{d}$ normalized eigenvector of $K^{-1}$:

(3.12)
$\left|\left|D_{opt}\right|\right|^{2}_{\mathcal{F}}=d^{2}\,\mathrm{Tr}\left(D_{opt}^{2}\right)=d^{2}||v||_{2}^{4}=d^{2}.$

For estimating $\left|\left|C_{opt}\right|\right|_{\mathcal{F}}$ we recall
$C_{opt}=D_{opt}K^{-1}+K^{1/2}\widetilde{J}_{opt}K^{-1/2}$, which implies
using the inequality
$||A\,B||_{\mathcal{F}}\leq||A||_{\mathcal{F}}||B||_{\mathcal{B}(\mathbb{R}^{d})}$
(see [10], p. 364):

(3.13)
$\left|\left|C_{opt}\right|\right|_{\mathcal{F}}\leq\left|\left|D_{opt}\right|\right|_{\mathcal{F}}\max(\sigma(K^{-1}))+\sqrt{\kappa(K)}\left|\left|\widetilde{J}_{opt}\right|\right|_{\mathcal{F}}\,.$

Since the Frobenius norm is unitarily invariant and
$\widetilde{J}_{opt}=\Psi\widehat{J}_{opt}\Psi^{-1}$ we have
$\left|\left|\widetilde{J}_{opt}\right|\right|_{\mathcal{F}}=\left|\left|\widehat{J}_{opt}\right|\right|_{\mathcal{F}}$.
For any $k=1,...,d$ we define $\alpha_{k}:=\langle v,\psi_{k}\rangle$ and we
observe that $\alpha_{k}^{2}=\frac{1}{d}$: Indeed from (3.5) it follows that

$\lambda_{opt}=\langle\psi_{k},\widetilde{D}_{opt}\psi_{k}\rangle=d\lambda_{opt}\alpha_{k}^{2}.$

Hence we can rewrite (3.6) as

(3.14)
$\left(\widehat{J}_{opt}\right)_{j,k}=\frac{\lambda_{j}+\lambda_{k}}{\lambda_{j}-\lambda_{k}}d\lambda_{opt}\alpha_{j}\alpha_{k}\,,\qquad\forall
j\neq k.$

It follows that

(3.15)
$\left|\left|\widehat{J}_{opt}\right|\right|_{\mathcal{F}}^{2}=d^{2}\lambda_{opt}^{2}\sum_{j\neq
k=1}^{d}\left(\frac{\lambda_{j}+\lambda_{k}}{\lambda_{j}-\lambda_{k}}\right)^{2}\alpha_{j}^{2}\alpha_{k}^{2}=\lambda_{opt}^{2}\sum_{j\neq
k=1}^{d}\left(\frac{\lambda_{j}+\lambda_{k}}{\lambda_{j}-\lambda_{k}}\right)^{2}.$

Next we choose the parameters $\lambda_{k}$, $k=1,...,d$ as

(3.16) $\lambda_{k}:=\frac{d-1}{c^{2}-1}+k-1,$

and they satisfy $0<\lambda_{1}<\cdots<\lambda_{d}$ and
$\frac{\lambda_{d}}{\lambda_{1}}=c^{2}$. Moreover we have for $j\neq k$:
$(\lambda_{j}+\lambda_{k})^{2}<(2\lambda_{d})^{2}=4\left(\frac{c^{2}}{c^{2}-1}\right)^{2}(d-1)^{2}$,
which implies with (3.15):

(3.17)
$\left|\left|\widehat{J}_{opt}\right|\right|_{\mathcal{F}}^{2}\leq\lambda_{opt}^{2}4\left(\frac{c^{2}}{c^{2}-1}\right)^{2}(d-1)^{2}\sum_{j\neq
k=1}^{d}\frac{1}{(\lambda_{j}-\lambda_{k})^{2}}.$

With the following estimate of a hyperharmonic series

$\sum_{j\neq
k=1}^{d}\frac{1}{(\lambda_{j}-\lambda_{k})^{2}}=\sum_{j=1}^{d}\sum_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{d}\frac{1}{(j-k)^{2}}=\sum_{j=1}^{d}\left(\sum_{l=1}^{j-1}\frac{1}{l^{2}}+\sum_{l=1}^{d-j}\frac{1}{l^{2}}\right)\leq\sum_{j=1}^{d}\frac{\pi^{2}}{3}=d\frac{\pi^{2}}{3},$

and (3.17) we obtain

(3.18)
$\left|\left|\widehat{J}_{opt}\right|\right|_{\mathcal{F}}\leq\lambda_{opt}\frac{2\pi}{\sqrt{3}}\frac{c^{2}}{c^{2}-1}\sqrt{d}(d-1).$

We conclude the proof by combining the (in)equalities (3.12), (3.13), and
(3.18). ∎

Let us briefly compare the strategy of proofs for Theorem 3.1 here and for
Theorem 2.2 in [8]: The main difference concerns how to connect the evolution
of the drift ODE to the FP-equation (here via the equality of the propagator
norms, and via a hypocoercive entropy method in [8]). Further, our choice of
the parameters $\lambda_{k}$ is (slightly) improved compared to the choice

(3.19) $\lambda_{k}=d+k,$

in [11, Remark 7] and [8]. Finally, the proof of Theorem 3.1(b) provides a
refined estimate of $\left|\left|C_{opt}\right|\right|_{\mathcal{F}}$.

######  Remark 3.2.

We note that, for any $c>1$, an optimal matrix pair
$(C_{opt}(c),\,D_{opt}(c))$ is not unique: Using in the proof of Theorem 3.1
the matrices $\widetilde{J}_{opt}^{T}$, $\widetilde{C}_{opt}^{T}$ instead of,
respectively, $\widetilde{J}_{opt}$, $\widetilde{C}_{opt}$ and the norm
$||\cdot||_{Q}$ instead of $||\cdot||_{P}$ yields another non-symmetric FP-
equation that satisfies the same estimates (3.1), (3.2).

##  4\. Examples and numerical illustrations

In this section we shall illustrate the results of §3. For an explicit example
in $\mathbb{R}^{2}$ we shall give a plot of the exact propagator norm for the
FP-equation, which is accessible due to Theorem 2.4 for constant-in-time
coefficients and due to Theorem 2.8 for the time-dependent case. First of all
we shall illustrate Theorem 3.1(a), particularly focussing on the
multiplicative constant in the exponential decay estimate (3.1).

###  4.1. Optimal decay estimates

As a first example we  consider the covariance matrix
$K=\mathrm{diag}(1,2)\in\mathbb{R}^{2\times 2}$. Then the maximum decay rate
for FP-equations that converge to $f_{\infty,K}$ is
$\lambda_{opt}=\min(\sigma(K))^{-1}=1$. Next we shall construct one optimal
pair of matrices $(C_{opt},D_{opt})$ such that
$e^{-L_{C_{opt},D_{opt}}t}f_{0}$ converges to $f_{\infty,K}$ with decay rate
$\lambda_{opt}$ and with a multiplicative constant arbitrary close to one. For
any $c>1$ we choose real numbers $0<\lambda_{1}<\lambda_{2}$ such that
$\frac{\lambda_{2}}{\lambda_{1}}=c^{2}$. We abbreviate
$\mu:=\frac{\lambda_{2}+\lambda_{1}}{\lambda_{2}-\lambda_{1}}=\frac{c^{2}+1}{c^{2}-1}>1$.
Following the procedure described in the proof of Theorem 3.1(a) we first
compute $D_{opt}=\widetilde{D}_{opt}=\mathrm{diag}(2,0)\in\mathbb{R}^{2\times
2}$. An orthonormal basis of $\mathbb{R}^{2}$ satisfying condition (3.5) is
given by $\psi_{1}:=\frac{1}{\sqrt{2}}(1,1)^{T}$ and
$\psi_{2}:=\frac{1}{\sqrt{2}}(-1,1)^{T}$. This defines the anti-symmetric
matrix  $\widetilde{J}_{opt}=\widehat{J}_{opt}=\begin{pmatrix}0&\mu\\\
-\mu&0\end{pmatrix}.$ Finally we compute

(4.1) $C_{opt}=\begin{pmatrix}2&\frac{\mu}{\sqrt{2}}\\\
-\sqrt{2}\mu&0\end{pmatrix},\text{ \quad and
\quad}\widetilde{C}_{opt}=\begin{pmatrix}2&\mu\\\ -\mu&0\end{pmatrix}.$

The spectral gaps of the drift matrix $C_{opt}$ and the operator
$L_{C_{opt},D_{opt}}$ coincide and are equal to 1. $C_{opt}$ has the two
distinct eigenvalues $1\pm i\sqrt{\mu^{2}-1}$ (because $\mu>1$), which are
also eigenvalues of $L_{C_{opt},D_{opt}}$ (see Theorem 5.3 in [5] or
Proposition 10 in [11]). Hence, $\lambda_{opt}=1$ is indeed the largest
possible, uniform decay rate of the FP-propagator $e^{-L_{C_{opt},D_{opt}}t}$
on $V_{0}^{\perp}$.

Thanks to Theorem 2.4 we can reduce the evaluation of the multiplicative
constant $c$ in the decay estimate (1.10) to the study of the propagator norm
of the associated drift ODE $\frac{d}{dt}{x}=-\widetilde{C}_{opt}x$. In
Theorem 3.7, [2] the authors provide the explicit form of the best
multiplicative constant for an ODE $\frac{d}{dt}{x}=-Ax$ in $\mathbb{R}^{2}$
when the matrix $A\in\mathbb{R}^{2\times 2}$ is positive stable,
diagonalizable and $\operatorname{Re}\tau_{1}=\operatorname{Re}\tau_{2}$, with
$\tau_{j}$, $j=1,2$ the eigenvalues of $A$: Then the best constant $c_{min}$
in the exponential decay estimate for $e^{-At}$ is given by

(4.2)
$c_{min}=\sqrt{\frac{1+\alpha}{1-\alpha}},\qquad\alpha:=\left|\left\langle\frac{v_{1}}{||v_{1}||_{2}},\frac{v_{2}}{||v_{2}||_{2}}\right\rangle\right|,$

where $v_{i}\in\mathbb{C}^{2}$, $i=1,2$ denote the eigenvectors of $A$. Since
$\widetilde{C}_{opt}$ satisfies the hypotheses of Theorem 3.7 in [2], a
straightforward computation gives $\alpha=\frac{1}{|\mu|}$, and the best
multiplicative constant is $c_{min}=\sqrt{\frac{\lambda_{2}}{\lambda_{1}}}=c$,
coinciding with the statement of Theorem 3.1(a).

We observe that $c\searrow 1$ implies $\mu\to\infty$. This corresponds to the
high-rotational limit in the drift matrix of the FP-equation

$\partial_{t}f_{t}=\mathrm{div}_{x}{(D_{opt}\nabla_{x}f_{t}+C_{opt}(\mu)xf_{t})}.$

For increasing $|\mu|$, the latter tends to mix with increasing speed the
dissipative and non-dissipative directions (i.e. $x_{1}$ and $x_{2}$,
respectively) of the corresponding symmetric FP-equation (i.e. with $\mu=0$).

As stated in Remark 3.2, replacing $\mu$ by $-\mu$ yields another FP-equation
with the same optimal decay behavior. Only the rotational direction is then
reversed.

###  4.2. Numerical illustrations: time-independent FP-equations

To illustrate the construction of optimal coefficient matrices in Theorem
3.1(a) we revisit the 2D-example from [8], i.e.
$K=\operatorname{diag}(1/\varepsilon,1)$, $\varepsilon=0.05$ which admits the
optimal decay rate $\lambda_{opt}=1$. For any given multiplicative constant
$c>1$, the optimal coefficient matrices constructed in Theorem 3.1(a) read:

(4.3) $D_{opt}=\widetilde{D}_{opt}=\begin{pmatrix}0&0\\\
0&2\end{pmatrix},\quad
C_{opt}=\begin{pmatrix}0&-\frac{\mu}{\sqrt{\varepsilon}}\\\
\sqrt{\varepsilon}\mu&2\end{pmatrix},\quad\widetilde{C}_{opt}=\begin{pmatrix}0&-\mu\\\
\mu&2\end{pmatrix},\quad\mu:=\frac{c^{2}+1}{c^{2}-1}.$

Figure 1. The solid curves show the FP- and ODE-propagator norms as functions
of $t$ for 3 values of the multiplicative parameter: $c=3,\,2,\,1.5$ (top to
bottom). The dashed curves give the corresponding (sharp) exponential bound of
the form $ce^{-\lambda_{opt}t}$ for the 3 cases. The dashed black curve shows
the exponential bound in the high-rotational limit, i.e. for $c\searrow 1$.
Colors only online.

In Figure1 we present the exact propagator norms (as a function of time) of
the FP-equation and of its drift ODE, i.e.

(4.4)
$\left\|e^{-L_{C,D}t}\right\|_{\mathcal{B}(V_{0}^{\perp})}=\left\|e^{-\widetilde{C}t}\right\|_{\mathcal{B}(\mathbb{R}^{d})},\quad
t\geq 0$

for several prescribed values of the multiplicative constant: $c=1.5,\,2,\,3$.
This figure includes also the r.h.s. of the corresponding exponential decay
estimate (3.1), using a logarithmic scale for the ordinate axis. Being the
exact upper envelops, this reveals that this estimate is indeed sharp,
concerning both the exponential rate and the multiplicative constant. Also
note that each curve of the propagator norm periodically touches (from above)
the curve corresponding to the high-rotational limit, given by
$e^{-\lambda_{opt}t}$.

Continuing with the same example, we shall next compare the results from
Theorem 3.1(a) here and Theorem 2.2 in [8]. First we need to explain the
criterion of comparison: For both results, and for a given constant $c>1$ we
seek a matrix pair $(C,D)$ such that the inequality (3.1) holds. Since
$c_{inf}=1$, such a pair can always be found, but $\|C\|_{\mathcal{F}}$
becomes large as $c\searrow 1$ (see §4.1). So, asking (4.4) to be close to the
high-rotational limit $e^{-\lambda_{opt}t}$ cannot be a useful criterion.
Instead, for given $c>1$ we want to find $(C,D)\in\mathcal{I}(K)$ such that
(3.1) holds and $\|C\|_{\mathcal{F}}$ is minimal. This has also a practical
implication for solving the FP-equation (2.2) numerically:
$\|C\|_{\mathcal{F}}$ “small” allows to use “large” time steps.
For fixed $c=\sqrt{2}$, Theorem 3.1(a) here and Theorem 2.2 in [8] yield,
respectively,

$\widetilde{C}_{opt}^{AS}=\begin{pmatrix}0&-3\\\
3&2\end{pmatrix},\quad\widetilde{C}_{opt}^{GM}=\begin{pmatrix}0&-7\\\
7&2\end{pmatrix},$

with $\|{C}_{opt}^{AS}\|_{\mathcal{F}}=\sqrt{184.45}$ and
$\|{C}_{opt}^{GM}\|_{\mathcal{F}}=\sqrt{986.45}$ . The essential difference
stems from the different choices of $\lambda_{1}$ and $\lambda_{2}$, (3.16)
vs. (3.19). In Theorem 3.1(a), the estimate (3.1) is sharp, and hence the
corresponding plot of the propagator norm has $\sqrt{2}\,e^{-t}$, i.e. the
r.h.s. in (3.1), as its upper envelop (see Figure 2, left). Since the estimate
from Theorem 2.2 in [8] is not sharp, the anti-symmetric part of
$\widetilde{C}_{opt}^{GM}$ is larger than “necessary”, and hence the
corresponding plot of the propagator norm stays well below the estimate
$\sqrt{2}\,e^{-t}$. With a view towards numerical applications the latter is
rather disadvantageous.

Figure 2. Left: For $c=\sqrt{2}$, the solid blue and red curves show the FP-
and ODE-propagator norms as functions of $t$ for the hypocoercive FP-equations
constructed, respectively, in Theorem 3.1(a) here and Theorem 2.2 in [8]. The
dashed blue curve gives the corresponding exponential bound
$\sqrt{2}\,e^{-t}$; it is sharp for Theorem 3.1(a). The solid green and black
curves show the FP- and ODE-propagator norms for the symmetric FP-equations in
(2.9) and (4.8), respectively.
Right: a zoom of the plot, close to $t=0$. Colors only online.

Figure 2 also shows the decay of the propagator norm of the symmetric FP-
equations in (2.9) and (4.8). Their respective decay rates are
$\frac{1}{\lambda_{max}(K)}=\varepsilon$ and
$\frac{d}{\operatorname{Tr}(K)}=\frac{2\varepsilon}{1+\varepsilon}$, both well
below $\lambda_{opt}=1$, the rate of the optimal hypocoercive FP-equations.

###  4.3. Numerical illustrations: time-dependent FP-equations

In [8] the authors used a FP-equation of the split form (2.9) with piecewise
constant coefficient matrices in order to approach the given equilibrium
quickly. Following this approach, we shall next discuss if time-dependent
coefficient matrices $C(t)$, $D(t)$ can accelerate the convergence in FP-
equations, compared to the case of constant matrices $C$, $D$ that was
analyzed in §3.
As a first step we shall show that the initial decay of hypocoercive FP-
evolutions, as constructed in the proof of Theorem 3.1(a) (recall that
$\operatorname{rank}(\widetilde{D}_{opt})=1$), can always be improved, e.g. in
the spirit of the split FP-equation (2.9) proposed in [8]. The following lemma
gives, at $t=0$, the largest possible decay rate of the FP-equation (2.2) as
well as of its drift ODE $\dot{x}=-\widetilde{C}x$ (both when considering
their propagator norms).

######  Lemma 4.1.

Let $K\in\mathcal{S}^{>0}$ be given. For any $(C,D)\in\mathcal{I}(K)$, the
maximum decay rate of
$\left\|e^{-L_{C,D}t}\right\|_{\mathcal{B}(V_{0}^{\perp})}$ at $t=0$ equals
$\frac{d}{\operatorname{Tr}(K)}$. It is obtained by the symmetric FP-equation
with $C=\frac{d}{\operatorname{Tr}(K)}\,I_{d}$ and
$D=\frac{d}{\operatorname{Tr}(K)}\,K$.

###### Proof.

Due to Theorem 2.4 we want to maximize the decay of the corresponding ODE-
propagator norm,

$\left\|e^{-\widetilde{C}t}\right\|_{\mathcal{B}(\mathbb{R}^{d})}^{2}=\lambda_{max}\left(e^{-\widetilde{C}^{T}t}e^{-\widetilde{C}t}\right)$

at $t=0$. A Taylor expansion yields

(4.5)
$\left\|e^{-\widetilde{C}t}\right\|_{\mathcal{B}(\mathbb{R}^{d})}=1-\lambda_{min}(\widetilde{C}_{s})\,t+\mathcal{O}(t^{2})\quad\mbox{
as }t\to 0,$

where $\widetilde{C}_{s}:=\frac{1}{2}(\widetilde{C}+\widetilde{C}^{T})$ is the
symmetric part of $\widetilde{C}$. We recall from the proof of Theorem 3.1(a)
that $\widetilde{C}:=K^{-1/2}CK^{1/2}$ and
$\widetilde{D}:=K^{-1/2}DK^{-1/2}=\widetilde{C}_{s}\geq 0$.
Thus we are led to the following optimization problem: Find
$\widetilde{C}_{s}\in\mathcal{S}^{\geq 0}$ with

(4.6)
$\operatorname{Tr}(D)=\operatorname{Tr}(K^{1/2}\widetilde{C}_{s}K^{1/2})=:\tau\leq
d,$

such that $\lambda_{min}(\widetilde{C}_{s})$ is maximal. Since
$\widetilde{J}$, the anti-symmetric part of $\widetilde{C}$, does not appear
within this problem, we set it to 0, for simplicity.
For such an optimal $\widetilde{C}_{s}$, (4.6) actually has to be an equality:
Otherwise we would have

$K^{1/2}\widetilde{C}_{s}K^{1/2}\leq\tau\,I_{d}<d\,I_{d}$

and the matrix $\widetilde{C}_{s}$ could be “enlarged”, e.g. by the matrix

$A:=\frac{d-\tau}{d^{2}-\tau}(d\,K^{-1}-\widetilde{C}_{s})\in\mathcal{S}^{>0}.$

Then, $\widetilde{C}_{s}+A$ still satisfies the constraint (4.6):

$\operatorname{Tr}\left(K^{1/2}[\widetilde{C}_{s}+A]K^{1/2}\right)=d,$

but $\lambda_{min}(\widetilde{C}_{s}+A)>\lambda_{min}(\widetilde{C}_{s})$,
contradicting the optimality of $\widetilde{C}_{s}$.
Next we shall prove that the optimal matrix satisfies

(4.7)
$\widetilde{C}_{s}=\widetilde{C}=\frac{d}{\operatorname{Tr}(K)}\,I_{d}=C.$

If the optimal $0\neq\widetilde{C}_{s}\in\mathcal{S}^{\geq 0}$ was not
proportional to $I_{d}$, we could “reduce” $\widetilde{C}_{s}$ by the matrix

$B:=\widetilde{C}_{s}-\lambda\,I_{d}\geq 0\quad\mbox{ with
}\lambda:=\lambda_{min}(\widetilde{C}_{s}),$

without changing the smallest eigenvalue. Moreover
$\widetilde{C}_{s}-B=\lambda\,I_{d}$ satisfies

$\operatorname{Tr}\left(K^{1/2}[\lambda\,I_{d}]K^{1/2}\right)=\operatorname{Tr}(K^{1/2}\widetilde{C}_{s}K^{1/2})-\operatorname{Tr}(K^{1/2}BK^{1/2})\leq
d,$

and hence $\lambda\,I_{d}$ is another optimal matrix of the above optimization
problem. From the equality requirement in (4.6) and
$\operatorname{Tr}(K^{1/2}\widetilde{C}_{s}K^{1/2})=d$ we then conclude $B=0$.
Hence $\widetilde{C}_{s}$ is proportional to $I_{d}$, and equality in (4.6)
yields $\lambda_{min}(\widetilde{C}_{s})=\frac{d}{\operatorname{Tr}(K)}$,
finishing the proof. ∎

With this lemma we can identify the symmetric FP-equation with steady state
$f_{\infty,K}$ that exhibits maximum initial decay as

(4.8)
$\partial_{t}f_{t}=\frac{d}{\operatorname{Tr}(K)}\mathrm{div}_{x}{(K\nabla_{x}f_{t}+xf_{t})},\qquad
x\in\mathbb{R}^{d},\ t\in(0,\infty).$

Its initial decay rate, $\frac{d}{\operatorname{Tr}(K)}$ is larger then that
of (2.9), namely $\frac{1}{\lambda_{max}(K)}$. We recall that the optimal FP-
equations constructed in the proof of Theorem 3.1(a) are all hypocoercive,
satisfying $\operatorname{rank}(\widetilde{D}_{opt})=1$, where
$\widetilde{D}_{opt}=(\widetilde{C}_{opt})_{s}$. Hence
$\lambda_{min}\big{(}(\widetilde{C}_{opt})_{s}\big{)}=0$, and the
corresponding propagator norm behaves like $1+\mathcal{O}(t^{2})$, see (4.5).
Therefore it is obvious that, for small time, the symmetric FP-equations (2.9)
and (4.8) both decrease the FP-propagator norm faster than the hypocoercive
FP-evolutions from Theorem 3.1(a). This is illustrated on a 2D example in
Figure 2, right.

For the rest of this section we shall base our discussion of using time-
dependent coefficients on the concrete example from §4.2, again with
$\varepsilon=0.05$, since a general theory of it seems unreachable to us for
the moment. In a numerical case study we shall analyze the FP-propagator norm
$\|S(t,0)\|_{\mathcal{B}(V_{0}^{\perp})}$, as a function of time. In the past
it would have been quite a challenge to compute (not just to estimate) this
norm. But due to Theorem 2.8 this has become easy for FP-equations with linear
drift.

Figure 3. Left: For $c=\sqrt{4/3}$, the FP- and ODE-propagator norms are given
for hypocoercive FP-equations with piecewise constant coefficients, using 5
different values on $0\leq t\leq 0.1$: The solid red curve corresponds to the
optimal, constant matrices from Theorem 3.1(a) as reference case, and the
dashed red curve is the corresponding decay estimate (3.1). The initially
symmetric FP-equations from (4.8) and (2.9) are given by the black and green
solid curves, respectively. Hypocoercive FP-equations with slower and faster
rotational drift are represented, respectively, by the blue and magenta solid
curves.
Right: a zoom of the plot, close to $t=0$. Colors only online.

In Figure 3 we shall compare the decay of the FP- and corresponding ODE-
propagator norms for 5 cases of FP-equations with piecewise constant
coefficient matrices, as in (2.9):

$S(t,0)=\begin{cases}e^{-L_{C_{i},D_{i}}t}&0\leq t\leq t_{0}\\\
e^{-L_{C_{1},D_{1}}(t-t_{0})}e^{-L_{C_{i},D_{i}}t_{0}}&t>t_{0}\end{cases},\quad
T(t,0)=\begin{cases}e^{-\widetilde{C}_{i}t}&0\leq t\leq t_{0}\\\
e^{-\widetilde{C}_{1}(t-t_{0})}e^{-\widetilde{C}_{i}t_{0}}&t>t_{0}\end{cases}.$

Choosing $t_{0}=0.1$, we use on the interval $(t_{0},\infty)$ always the same
matrices, namely those from (4.3) with $\mu=7$, which is the optimal
hypocoercive FP-evolution from Theorem 3.1 for the multiplicative constant
$c=\sqrt{4/3}$. For the interval $[0,t_{0}]$ we shall compare the following
cases:

  1. (FP1)

This reference case uses the same coefficients as for $t>t_{0}$, i.e.:

(4.9) $D_{1}=\operatorname{diag}(0,2),\quad
C_{1}=\begin{pmatrix}0&-\frac{7}{\sqrt{\varepsilon}}\\\
7\sqrt{\varepsilon}&2\end{pmatrix},\quad\widetilde{C}_{1}=\begin{pmatrix}0&-7\\\
7&2\end{pmatrix}.$

Figure 3 also includes the sharp upper and lower envelops of the resulting
propagator norm (as function of $t$).

  2. (FP2)

The symmetric FP-equation from [8], and shown in (2.9) is determined by the
matrices

$D_{2}=I_{2},\quad
C_{2}=\widetilde{C}_{2}=\operatorname{diag}(\varepsilon,1).$

  3. (FP3)

The symmetric FP-equation (4.8) with maximum initial decay is determined by
the matrices

$D_{3}=\frac{2\varepsilon}{1+\varepsilon}\operatorname{diag}(\frac{1}{\varepsilon},1),\quad
C_{3}=\widetilde{C}_{3}=\frac{2\varepsilon}{1+\varepsilon}I_{2}.$

  4. (FP4)

A hypocoercive FP-equation with slower rotational part than in (4.9) is
determined by the matrices

$D_{4}=\operatorname{diag}(0,2),\quad
C_{4}=\begin{pmatrix}0&-\frac{3}{\sqrt{\varepsilon}}\\\
3\sqrt{\varepsilon}&2\end{pmatrix},\quad\widetilde{C}_{4}=\begin{pmatrix}0&-3\\\
3&2\end{pmatrix}.$

  5. (FP5)

A hypocoercive FP-equation with faster rotational part than in (4.9) is
determined by the matrices

$D_{5}=\operatorname{diag}(0,2),\quad
C_{5}=\begin{pmatrix}0&-\frac{11}{\sqrt{\varepsilon}}\\\
11\sqrt{\varepsilon}&2\end{pmatrix},\quad\widetilde{C}_{5}=\begin{pmatrix}0&-11\\\
11&2\end{pmatrix}.$

Note that (FP4) and (FP5) are both of the form (4.3).

First we need to fix the criterion for comparing these 5 FP-equations with
split coefficients. As one sees from Figure 3, adapting the FP-equation only
on the initial time interval $[0,t_{0}]$ has a highly nonlocal-in-$t$ effect.
Hence, it does not make sense to compare the norm-curves pointwise in time.
Following the paradigm of §3, it is appropriate to compare again the
corresponding best exponential decay estimates (3.1). Since all compared FP-
equations coincide for large time, or more precisely on $(t_{0},\infty)$,
their exponential decay rate is the same, and it suffices to compare the
multiplicative constant of the (sharp) decay estimates.
Now we shall replace in the reference FP-equation (FP1) the initial phase by a
symmetric evolution: With both options (FP2) and (FP3) the propagator norm
decays initially faster than for the reference FP-equation (see Figure 3,
right), but this backfires at later times: In both cases the upper envelop for
the whole norm-function on $[0,\infty)$ and hence the multiplicative constant
$c$ is larger than for the reference case (FP1) (see Figure 3, left).
Finally we shall replace in the reference FP-equation (FP1) the initial phase
by a hypocoercive evolution having an anti-symmetric part of $\widetilde{C}$
that differs from case (FP1). With the slower rotational part in case (FP4)
the multiplicative constant $c$ is increased (see Figure 3, left), but when
using initially the faster rotational part from case (FP5), the multiplicative
constant $c$ is decreased!

Figure 4. For $c=\sqrt{4/3}$, the FP- and ODE-propagator norms are given for
hypocoercive (h.c.) FP-equations with piecewise constant coefficients, using 3
different values on $0\leq t\leq t_{0}$: The solid red curve corresponds to
the optimal, constant matrices from Theorem 3.1(a) as reference case.
Hypocoercive FP-equations with the faster rotational drift matrices (FP5),
(FP6) are represented by the magenta and blue solid curves, respectively. The
dashed curves are the corresponding decay estimate (3.1). The discontinuity
points $t_{0}$ of the coefficient matrices are marked with black dots.

While we present in Figure 3 the plots only for $t_{0}=0.1$, the results for
other values of $t_{0}>0$ are qualitatively the same. Choosing $t_{0}\approx
0.1434$ (i.e. the first point of tangency between
$\|e^{-\widetilde{C}_{5}t}\|_{\mathcal{B}(\mathbb{R}^{d})}$ and its sharp
exponential decay estimate $\sqrt{6/5}e^{-t}$, see Figure 4) in the split case
(FP5) reduces the multiplicative constant to $c=\sqrt{6/5}$. Note that this is
also the sharp constant for the non-split FP-equation involving the matrices
$(C_{5},\,D_{5})$. This means that the same decay quality (in the above
defined sense) can be obtained with the constant coefficient matrices
$(C_{5},\,D_{5})$ for all time or just a short initial layer with
$(C_{5},\,D_{5})$ and then evolving with $(C_{1},\,D_{1})$ for $t>t_{0}$. The
multiplicative constant can be reduced even further, e.g. with the following
choice of matrices on the interval $[0,0.11413]$ (see Figure 4):

  1. (FP6)

$D_{6}=\operatorname{diag}(0,2),\quad
C_{6}=\begin{pmatrix}0&-\frac{13.8}{\sqrt{\varepsilon}}\\\
13.8\sqrt{\varepsilon}&2\end{pmatrix},\quad\widetilde{C}_{6}=\begin{pmatrix}0&-13.8\\\
13.8&2\end{pmatrix}$

This example of time-dependent FP-coefficients is also algorithmically
relevant, since $\|C_{1}\|_{\mathcal{F}}<\|C_{6}\|_{\mathcal{F}}$. Hence,
longer time steps could be used in a discretization of the split FP-equation
for $t>t_{0}$.

##  5\. Conclusion

For any given anisotropic Gaussian steady state (2.1) with covariance matrix
$K$, we analyzed the construction of non-symmetric FP-equations (2.2) that
show fastest decay towards the unique normalized steady state $f_{\infty,K}$.
Building upon preceding results (in particular [11, 8]) we proved that optimal
exponential decay with small multiplicative constants (as in (1.10), and
uniformly in $f_{0}$) can be achieved with a single FP-equation, without
having to split off an initial evolution phase. Thereby, the maximum decay
rate $\lambda_{opt}=\max(\sigma(K^{-1}))$, and the infimum of the
multiplicative constants $c_{inf}=1$. By contrast, the best multiplicative
constant obtainable in [8] was bounded below by $\sqrt{\kappa(K)\,e}$. Hence,
the gain provided here for the multiplicative constant is particularly
important when $\kappa(K)$ is large, i.e. when the original, symmetric FP-
dynamics includes very different time scales due to very different eigenvalues
in $K$.

More precisely, for any given multiplicative constant $c>1$ we were able to
construct explicitly a non-symmetric FP-equation of form (2.2) with constant
drift matrix $C_{opt}(c)$ and diffusion matrix $D_{opt}(c)$ such that the
exponential decay estimate (1.10) holds with the parameters
$(\lambda_{opt},\,c)$. For given $c$ and variable space dimension $d$, we were
able to reduce the growth estimate on these drift matrices to
$\mathcal{O}(d^{3/2})$, down from $\mathcal{O}(d^{2})$ given in [8].

In explicit 2D examples we illustrated, both analytically and numerically,
that the infimum of the multiplicative constant, $c_{inf}=1$ corresponds to
the limit of adding a highly rotational, non-symmetric drift to the original
FP-equation.

To round off our analysis we presented a numerical case study on a FP-equation
in 2D with piecewise constant coefficient matrices. This showed two unexpected
phenomena: First, no symmetric FP-evolution on an initial time layer was able
to improve the overall decay behavior; in fact it always got worse than in the
time-independent case. Second, replacing on an initial time layer the non-
symmetric drift by a higher rotational one (and then returning to the original
drift for all time) can reduce the multiplicative constant for the whole
evolution to a level that pertains to a “larger” drift matrix $C$.

##  Appendix A Proof of Theorem 2.8

###### Proof-idea.

First, the coordinate transformation $\tilde{x}:=K^{-1/2}x$ and
$\tilde{f}(\tilde{x}):=(\operatorname{det}K)^{1/2}\,f(K^{1/2}\tilde{x})$
transforms (2.12) into the normalized FP-equation

(A.1)
$\partial_{t}\tilde{f}_{t}=-\widetilde{L}(t)\tilde{f}_{t}:=\mathrm{div}_{\tilde{x}}{(\widetilde{D}(t)\nabla_{\tilde{x}}\tilde{f}_{t}+\widetilde{C}(t)\tilde{x}\tilde{f}_{t})},\qquad\tilde{x}\in\mathbb{R}^{d},\
t\in(0,\infty),$

where $\widetilde{D}(t):=K^{-1/2}D(t)K^{-1/2}$. This FP-equation is naturally
considered in
$\widetilde{\mathcal{H}}:=L^{2}(\mathbb{R}^{d},\tilde{f}_{\infty}^{-1})$, and
the (transformed) steady state is

$\tilde{f}_{\infty}(\tilde{x})=(2\pi)^{-d/2}\,e^{-|\tilde{x}|^{2}/2}.$

This transformation preserves the norm of the solution:
$\|f_{t}\|_{\mathcal{H}}=\|\tilde{f}_{t}\|_{\widetilde{\mathcal{H}}},\>t\geq
0$. Hence the propagator norms of (2.12) and (A.1) coincide:
$\left\|S(t_{2},t_{1})\right\|_{\mathcal{B}(V_{0}^{\perp})}=\left\|\widetilde{S}(t_{2},t_{1})\right\|_{\mathcal{B}(\widetilde{V}_{0}^{\perp})}$.
Next, one decomposes $\widetilde{\mathcal{H}}$ into mutually orthogonal
subspaces $\widetilde{V}^{(m)},\>m\in\mathbb{N}_{0}$, which are each invariant
under the operators $\widetilde{L}(t)$ $\forall t\geq 0$:

$\widetilde{\mathcal{H}}=\bigoplus_{m\in\mathbb{N}_{0}}{\\!\\!\\!}^{\perp}\,\widetilde{V}^{(m)},$

with

$\widetilde{V}^{(m)}:=\operatorname{span}\\{g_{\alpha}(\tilde{x}):=(-1)^{|\alpha|}\nabla^{\alpha}\tilde{f}_{\infty}(\tilde{x})\,:\,\alpha\in\mathbb{N}_{0}^{d},\,|\alpha|=m\\}.$

Decomposing the solution of (A.1) into these subspaces as

$\tilde{f}_{t}(\tilde{x})=\sum_{\alpha\in\mathbb{N}_{0}^{d}}\tilde{d}_{\alpha}(t)\,\frac{g_{\alpha}(\tilde{x})}{\|g_{\alpha}\|_{\widetilde{\mathcal{H}}}},$

yields the estimates

$\sum_{|\alpha|=m}|\tilde{d}_{\alpha}(t_{2})|^{2}\leq
h(t_{2},t_{1})^{2m}\left(\sum_{|\alpha|=m}|\tilde{d}_{\alpha}(t_{1})|^{2}\right),\quad
0\leq t_{1}\leq t_{2}<\infty,\>m\in\mathbb{N},$

with

$h(t_{2},t_{1}):=\|T(t_{2},t_{1})\|_{\mathcal{B}(\mathbb{R}^{d})}\leq 1,\quad
0\leq t_{1}\leq t_{2}<\infty.$

On the one hand this shows that

$\|\tilde{f}_{t_{2}}-\tilde{f}_{\infty}\|_{\widetilde{\mathcal{H}}}=\|\widetilde{S}(t_{2},t_{1})(\tilde{f}_{t_{1}}-\tilde{f}_{\infty})\|_{\widetilde{\mathcal{H}}}\leq\|T(t_{2},t_{1})\|_{\mathcal{B}(\mathbb{R}^{d})}\|\tilde{f}_{t_{1}}-\tilde{f}_{\infty}\|_{\widetilde{\mathcal{H}}},\quad
0\leq t_{1}\leq t_{2}<\infty.$

On the other hand we can use initial conditions
$\tilde{f}_{t_{1}}\in\widetilde{V}^{(1)}$, noting as in [4, §4.2] that the
coefficient vector
$\tilde{d}^{(1)}(t):=\left(\tilde{d}_{\alpha}(t)\right)_{|\alpha|=1}\in\mathbb{R}^{d}$
evolves according to

$\frac{d}{dt}\tilde{d}^{(1)}=-\widetilde{C}(t)\,\tilde{d}^{(1)},$

i.e. the drift ODE of the FP-equation. This implies the reverse inequality

$\left\|\widetilde{S}(t_{2},t_{1})\right\|_{\mathcal{B}(\widetilde{V}_{0}^{\perp})}\geq\left\|T(t_{2},t_{1})\right\|_{\mathcal{B}(\mathbb{R}^{d})},\qquad\forall
0\leq t_{1}\leq t_{2}<\infty,$

and hence the equality (2.14) follows. ∎

## Acknowledgement

The authors were partially supported by the FWF (Austrian Science Fund) funded
SFB #F65 and the FWF-doctoral school W 1245. We acknowledge fruitful
discussions with José A. Carrillo, who originally proposed this topic to us.
We are also grateful to the anonymous referees, whose suggestions helped to
improve this work.

## References

  * [1] F. Achleitner, A. Arnold, E. A. Carlen, On multi-dimensional hypocoercive BGK models. Kinetic and Related Models 11, No. 4 (2018), 953-1009.
  * [2] F. Achleitner, A. Arnold, B. Signorello, On optimal decay estimates for ODEs and PDEs with modal decomposition. Stochastic dynamics out of equilibrium, Springer Proceedings in Mathematics and Statistics 282 (2019), 241-264, G. Giacomin et al. (eds.), Springer.
  * [3] A. Arnold, P. A. Markowich, G. Toscani, A. Unterreiter, On convex Sobolev inequalities and the rate of convergence to equilibrium for Fokker-Planck type equations. Comm. PDE 26, No. 1-2 (2001), 43-100.
  * [4] A. Arnold, C. Schmeiser, B. Signorello, Propagator norm and sharp decay estimates for Fokker-Planck equations with linear drift. Preprint. https://arxiv.org/abs/2003.01405.
  * [5] A. Arnold, J. Erb, Sharp entropy decay for hypocoercive and non-symmetric Fokker-Planck equations with linear drift. Preprint. https://arxiv.org/abs/1409.5425.
  * [6] P. Diaconis, The Markov chain Monte Carlo revolution. Bulletin of the American Mathematical Society 46, no.2 (2009), 179-205.
  * [7] H. Dietert, J. Evans, Finding the jump rate for fastest decay in the Goldstein-Taylor model. Preprint. https://arxiv.org/abs/2103.10064.
  * [8] A. Guillin, P. Monmarché, Optimal linear drift for the speed of convergence of an hypoelliptic diffusion. Electron. Commun. Probab. 21, no. 74 (2016), 1-14; Erratum: Electron. Commun. Probab. 22, no. 15 (2017), 1-2.
  * [9] L. Miclo, P. Monmarché, Étude spectrale minutieuse de processus moins indécis que les autres, In: C. Donati-Martin, A. Lejay, A. Rouault (eds) Séminaire de Probabilités XLV, Lecture Notes in Mathematics, vol 2078, pp. 459–481. Springer, Heidelberg (2013). English summary available at https://www.ljll.math.upmc.fr/~monmarche
  * [10] R.A. Horn, C.R. Johnson, Matrix Analysis, 2nd edition, Cambridge University Press (2013).
  * [11] T. Lelièvre, F. Nier, G.A. Pavliotis, Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion. Journal of Statistical Physics, 152, no. 2 (2013), 237-274.
  * [12] J. Snyders, M. Zakai, On nonnegative solutions of the equation $AD+DA^{\prime}=-C$. SIAM J. Appl. Math. 18 (1970), 704-715.
  * [13] C. Villani, Hypocoercivity, Memoirs of the American Mathematical Society 202 (2009).


This paper proposes a parametric-based network architecture for joint channel estimation and data detection in communications systems with hardware impairments. This architecture is composed of a data-augmented layer, a custom soft thresholding function, and several linear layers modeling the effect of channel effects and hardware impairments. In the proposed network, the soft thresholding function softly constrains the detected data to be within the considered constellation. The latter depends only on one one parameter that is optimized during training. The benefit of the proposed approach is illustrated through a communication chain corrupted by multiple impairments and noises.

Communications Systems, Parametric Estimation, Machine Learning.

§ INTRODUCTION

In communications systems, the use of advanced coding schemes and modulation formats makes the performance of the system quite sensible to signal distortion occurring at the physical layer. These distortions include the effect of the propagation channel, non-ideal synchronization, IQ imbalance, and carrier impairments. In the literature, a large number of digital algorithms have been developed for hardware impairment compensation and data detection <cit.>. These algorithms are mainly based on parametric signal processing approaches. Recently, several studies have investigated the use of Deep Learning (DL) techniques to address some physical-layer problems <cit.>. In particular, many studies have focused on the use of model-based DL techniques for Multiple-Input Multiple-Output (MIMO) detection <cit.>. The main drawback of these approaches is that they usually require perfect channel state information. To address this issue, a DL network addressing the joint channel estimation and data detection problem has been proposed in <cit.>. As this technique does not assume any apriori model for the channel matrix, it often requires a large training database and is not suited for time-varying channels.

In this paper, we focus on the joint channel estimation and data detection problem in frequency-selective Single-Input Single-Output (SISO) communications. By modeling the channel as a cascade of parametric models, we propose a new parametric multi-layer network for channel estimation and data detection which is well suited for time-varying channels. This network can be trained using a small number of pilot symbols. Following the principle of Iterative thresholding algorithms, we propose to include a non-linear soft thresholding layer in the structure of the parametric network. This layer is based on the optimal denoiser for Gaussian noise that has been introduced in <cit.>.

The paper is organized as follows. Section <ref> describes the problem statement, Section <ref> presents the proposed network architecture, and Section <ref> illustrates the benefit of our approach together with the obtained simulation results.

[right, text centered,text width=4em, ] at (-4.25,1) $\mathbf{s}_p$
[->,>=latex] (-3.25,1)–(-2.5,1);
(-2.5,0.25) rectangle (-1,1.75) node[pos=.5,text width=4em, text centered](l1) $\mathbf{f}(.;\mathbf{u})$ ;
[->,>=latex] (-1,1)–(-0.25,1) node[midway,yshift=1em]$\mathbf{s}$ ;
(-0.25,0.25) rectangle (1.25,1.75) node[pos=.5,text width=4em, text centered](l1) $\mathcal{S}_{\sigma_s}(.)$ ;
[->,>=latex] (1.25,1)–(2,1) node[midway,yshift=1em]$\mathbf{s}_1$ ;
(2,0.25) rectangle (3.5,1.75) node[pos=.5,text width=4em, text centered](l1) $\mathbf{H}_1(\boldsymbol \theta_1)$ ;
[->,>=latex] (3.5,1)–(4.25,1) node[midway,yshift=1em]$\mathbf{s}_2$ ;
[dashed] (4.25,0.25) rectangle (5.75,1.75) node[pos=.5,text width=4em, text centered] $\cdots$;
[->,>=latex] (5.75,1)–(6.5,1);
(6.5,0.25) rectangle (8,1.75) node[pos=.5,text width=4em, text centered]  (l3)$\mathbf{H}_L(\boldsymbol \theta_L)$;
[->,>=latex] (8,1)–(8.75,1) node[right] $\widehat{\mathbf{y}}=\mathbf{H}(\boldsymbol \theta)\mathcal{S}_{\sigma_s}(\mathbf{f}(\mathbf{s}_p;\mathbf{u}))$;
at (-1.75,2.5) Data Augmented;
at (-1.75,2.2) Layer;
at (0.5,2.5) Non-linear Activation;
at (0.5,2.2) Function;
[right] at (-4.15,-0.25) learnable;
[right] at (-4.25,-0.75) parameters;
[right] at (-2.,-0.5) $\mathbf{u}$;
[right] at (0,-0.5) $\sigma_s^2$;
[right] at (2.5,-0.5) $\boldsymbol \theta_1$;
[right] at (4.75,-0.5)  $ \cdots$ ;
[right] at (7,-0.5) $\boldsymbol \theta_L$;

Proposed ParamNet network architecture. The network is composed of a data augmented layer, $\mathbf{f}(.;\mathbf{u})$, a non-linear activation function, $\mathcal{S}_{\sigma_s}(.)$, and $L$ linear layers $\mathbf{H}_l(\boldsymbol \theta_l)$ modeling the hardware impairments and channel effects.

§ PROBLEM STATEMENT

§.§ Signal Model

Let us consider a communications channel with $N$ complex-valued transmitted
symbols $\mathbf{s}_c=[s_{c}[0],\cdots,s_{c}[N-1]]^T$ and $N$ complex-valued received samples $\mathbf{y}_c=[y_{c}[0],\cdots,y_{c}[N-1]]^T$, where $^T$ corresponds to the vector transpose. Let us denote the augmented real-valued transmitted symbols and received samples by $\mathbf{s}=[\Re e(\mathbf{s}_c^T),\Im m(\mathbf{s}_c^T)]^T$ and $\mathbf{y}=[\Re e(\mathbf{y}_c^T),\Im m(\mathbf{y}_c^T)]^T$, respectively. By assuming a linear channel, the augmented received samples can be expressed by the following linear model:

\begin{align}
\mathbf{y} = \mathbf{H}(\boldsymbol \theta)\mathbf{s}+\mathbf{b},\label{eq1}
\end{align}

where $\mathbf{H}(\boldsymbol \theta)$ is a $2N\times 2N$ global transfer matrix, $\boldsymbol \theta$ corresponds to the unknown channel parameters, and $\mathbf{b}\sim \mathcal{N}(\mathbf{0},\sigma^2_b \mathbf{I})$ is an additive noise component with zero mean and variance $\sigma^2_b$. In communications systems, the transmitted data usually belong to a finite alphabet. This property implies that each augmented symbol $s[n]$ belongs to a set $\mathcal{M}=\{M_1,\cdots,M_K\}$ composed of $K$ real-valued symbols. Finally, we consider in this study that the global transfer matrix can be decomposed as a multi-layer channel as follows

\begin{align}
\mathbf{H}(\boldsymbol \theta)=\mathbf{H}_L(\boldsymbol \theta_L)\times \cdots \times \mathbf{H}_1(\boldsymbol \theta_1),
\end{align}

where $\mathbf{H}_l(\boldsymbol \theta_l)$ corresponds to the transfer matrix of the $l^{th}$ layer and depends on the unknown parameters $\boldsymbol \theta_l$. Note that this multi-layer model arises naturally in communications systems corrupted by multiple linear impairments or effects such as transmitter or receiver IQ imbalance, Carrier Frequency Offset, Phase Noise, frequency selective channels, etc.

§.§ Joint Estimation and Detection Problem

The objective of this paper is to jointly estimate $\boldsymbol \theta$ and $\mathbf{s}\in \mathcal{M}^{2N}$ from $\mathbf{y}$.
For Gaussian noise, this optimization problem can be formulated as $\min_{\boldsymbol \theta,\mathbf{s}\in \mathcal{M}^{2N}}\|\mathbf{y}-\mathbf{H}(\boldsymbol \theta)\mathbf{s}\|^2_2
$. This problem is challenging since $\boldsymbol \theta$ is unknown and $\mathbf{s}$ belongs to a discrete alphabet. To simplify this problem, one solution is to consider that $N_x$ complex-valued pilot symbols, $\mathbf{s}_{c}$, are known at the receiver. By using the knowledge of the real and imaginary parts of the pilot symbols, denoted by the augmented vector $\mathbf{s}_p=[\Re e(\mathbf{s}_{c}^T),\Im m(\mathbf{s}_{c}^T)]^T$, the optimization problem can be reformulated as

\begin{align}
\min_{\boldsymbol \theta,\mathbf{u}\in \mathcal{M}^{2N_u}}\|\mathbf{y}-\mathbf{H}(\boldsymbol \theta)\mathbf{f}(\mathbf{s}_p;\mathbf{u})\|^2_2,
\end{align}

where $\mathbf{u}$ is an unknown vector composed of $2N_u$ symbols, and $\mathbf{s}=\mathbf{f}(\mathbf{s}_p;\mathbf{u})$ is a data-augmented function that takes $2N_x$ real-valued symbols and returns $2N=2(N_x+N_u)$ real-valued symbols.

To estimate $\boldsymbol \theta$ and $\mathbf{u}$, two approaches can be employed. In the first approach, the alphabet of $\mathbf{u}$ can be simply neglected by assuming $\mathbf{u}\in \mathbb{R}^{2N_u}$. In the second approach, the alphabet can be enforced using a projected gradient algorithm by mapping periodically $\mathbf{u}$ to the closest point in $\mathcal{M}^{2N_u}$. Even if this second approach has the distinct advantage of exploiting the data constellation, the influence of projection errors in the medium-SNR region cannot be completely excluded. Furthermore, during training, the use of a hard-thresholding function makes it impossible to backpropagate the gradient to the data augmented layer since the function derivative is $0$ almost everywhere. To address these two issues, we propose a new network architecture composed of a soft projection layer.

§ PROPOSED MULTI-LAYER NETWORK

The structure of the proposed Multi-Layer parametric network is presented in Fig. <ref>. This section describes the mathematical models of each layer and presents the associated learning algorithm.

§.§ Layer Models

In this section, the real-valued inputs and outputs of each layer are denoted $\mathbf{x}_i$ and $\mathbf{x}_o$, respectively.

§.§.§ Data Augmented Layer

[scale=0.85, every node/.style=scale=0.85]
text width=2em,
minimum height=2em,
text centered,
[cell,fill=black!10] (l1) $x_i[0]$ ;
[cell,fill=black!10,below of=l1,node distance=2em] (l2) $x_i[1]$ ;
[cell,fill=black!10,below of=l2,node distance=2em] (l3) $x_i[2]$ ;
[cell,fill=black!10,below of=l3,node distance=2em] (l4) $\cdots$ ;
[below of=l4,node distance=3em,text width=6em] (l5) pilots $\mathbf{x}_i=\mathbf{s}_p$
[rectangle, draw, right of=l2,text width=5em, text centered,yshift=-1em,node distance=8em,minimum height=6em] (r1) Data Augmented Layer ;
[cell,fill=black!10,right of=l1,node distance=16em,yshift=2em] (x1) $x_i[0]$ ;
[cell,below of=x1,node distance=2em] (x2) $u[0]$ ;
[cell,below of=x2,node distance=2em] (x3) $u[1]$ ;
[cell,below of=x3,node distance=2em] (x4) $\cdots$ ;
[cell,fill=black!10,below of=x4,node distance=2em] (x5) $x_i[1]$ ;
[cell,below of=x5,node distance=2em] (x6) $\cdots$ ;
[below of=x6,node distance=3em,text width=6em] (l5) pilots + data
[left of=r1,node distance=6em] (r0);
[right of=r1,node distance=6em] (r2);
[->,>=latex] (r0)–(r1);
[->,>=latex] (r1)–(r2);

Illustration of the Data Augmented Layer. The learnable parameters correspond to the data  $\mathbf{u}=[u[0],u[1],\cdots]^T$.

The purpose of the data augmented layer is to insert some data $\mathbf{u}$ into the input vector $\mathbf{x}_i$ (see Fig. <ref>). Mathematically, the input-output relationship of this layer can be expressed by

\begin{align}
\mathbf{x}_o=\mathbf{f}(\mathbf{x}_i;\mathbf{u})=\mathbf{P}_0\mathbf{x}_i+\mathbf{P}_1\mathbf{u},
\end{align}

where the layer trainable parameters correspond to the unknown real-valued vector $\mathbf{u}$. The matrices $\mathbf{P}_0$ and $\mathbf{P}_1$ correspond to allocation matrices of size $2N\times 2N_x$ and $2N\times 2N_u$, respectively. An allocation matrix is a matrix with entries from $\{0,1\}$ that only contains a single one in each row.

§.§.§ Non-linear Activation Function

The non-linear activation layer is based on the soft thresholding operator

\begin{align}
\mathbf{x}_o=\mathcal{S}_{\sigma_s}(\mathbf{x}_i).
\end{align}

The purpose of the function $\mathcal{S}_{\sigma_s}(.)$ is to mimic a soft projection into the data constellation. To this end, we propose to use the optimal denoiser for Gaussian noise, which has been recently used in MIMO detection problems <cit.>. By assuming  an i.i.d. Gaussian distribution with diagonal covariance matrix $(\sigma_s^2/2) \mathbf{I}$ for the error $\mathbf{x}_i-\mathbf{s}$ and a uniform distribution over $\mathcal{M}$ for the data symbols $\mathbf{s}$, the $n^{th}$ output of the optimal denoiser can be expressed as
\begin{align}
[\mathcal{S}_{\sigma_s}(\mathbf{x}_i)]_n = \frac{ \sum_{s\in \mathcal{M}} s e^{-\frac{1}{\sigma_s^2}(x_i[n]-s)^2}}{\sum_{s\in \mathcal{M}} e^{-\frac{1}{\sigma_s^2}(x_i[n]-s)^2}}
\end{align}

This layer contains only one learnable parameter called the data noise variance $\sigma_s^2$.

For illustration purpose, Fig. <ref> depicts the function $\mathcal{S}_{\sigma_s}(.)$ for different values of $\sigma_s^2$. Note that the data noise variance $\sigma_s^2$ is different from the additive noise variance $\sigma_b^2$ in (<ref>). The data noise variance $\sigma_s^2$ allows to alternate between the exploration and exploitation of the data vector $\mathbf{u}$ during training. For example, the exploration is encouraged when $\sigma_s^2= 2$ since  $\mathcal{S}_{\sigma_s}(.)$ is nearly linear between $M_1$ and $M_4$, while the exploitation is encouraged when $\sigma_s^2=0.3$ since this function is close to a threshold detector.

[declare function =
S(,) = (-3*exp(-((x+3)^2)/)-1*exp(-((x+1)^2)/)+1*exp(-((x-1)^2)/)+3*exp(-((x-3)^2)/))/(exp(-((x+3)^2)/)+exp(-((x+1)^2)/)+exp(-((x-1)^2)/)+exp(-((x-3)^2)/));

[xlabel=$z$,ylabel=$\mathcal{S}_{\sigma_s}(z)$,xmin=-5,xmax=5,grid=both,width=0.9,legend pos=south east,xtick=-3,-1,1,3,xticklabels=$M_1$,$M_2$,$M_3$,$M_4$,
[black, thick,samples=100]S(x,0.01);
[blue, thick,samples=100] S(x,2);
[red, thick,samples=100]S(x,0.3);
Hard Projector, $\sigma_s^2=2$, $\sigma_s^2=0.3$

Soft projection into $\mathcal{M}=\{M_1,M_2,M_3,M_4\}$.

§.§.§ Parametric SISO Channel Layers

This section presents a list of linear layers commonly encountered in communications systems. While this list is not exhaustive and only focuses on SISO layers, it can be easily adapted to more specific problems or MIMO systems. For SISO layers, the augmented outputs can be expressed as $\mathbf{x}_{o}=\mathbf{H}_l(\boldsymbol \theta_l)\mathbf{x}_{i}$, where the augmented inputs $\mathbf{x}_{i}=[\Re e(\mathbf{x}_{c,i}^T),\Im m(\mathbf{x}_{c,i}^T)]^T$ contains the real and imaginary parts of the complex-valued input vector $\mathbf{x}_{c,i}=[x_{c,i}[0],\cdots,x_{c,i}[N-1]]^T$. Herein, the mathematical models are expressed in scalar form for sake of simplicity.

* IQ Imbalance Layer: The output of an IQ imbalance layer is given by

\begin{align}
x_{c,o}[n]=\mu x_{c,i}[n]+\nu x_{c,i}^*[n],
\end{align}

where $(.)^*$ corresponds to the complex conjugate, and $(\mu,\nu)\in\mathbb{C}^2$ are the IQ imbalance parameters <cit.>. This layer depends on the 4 real-valued parameters $\boldsymbol \theta_l=[\Re e(\mu),\Re e(\nu),\Im m(\mu), \Im m(\nu)]^T$.

* FIR Layer: The output of a FIR layer is given by

\begin{align}
x_{c,o}[n]=\sum_{k=0}^{D-1}h[k] x_{c,i}[n-k].
\end{align}

This layer depends on the real-valued vector $\boldsymbol \theta_l=[\Re e(\mathbf{h}^T),\Im m(\mathbf{h}^T)]^T$ where $\mathbf{h}=[h[0],\cdots,h[D-1]]^T$ corresponds to the channel impulse response.

* Piecewise-Constant Phase Layer: The output of this layer is given by

\begin{align}
x_{c,o}[n]=e^{j\varphi[n~\text{mod} N_s]} x_{c,i}[n]
\end{align}

where mod corresponds to the modulo operator <cit.>. This layer depends on $N/N_s$ real-valued parameters that are given by $\boldsymbol \theta_l=[\varphi[0],\cdots,\varphi[L-1]]^T$.

§.§ Network training/operating

To train the network, we propose to minimize the custom regularized loss function given by

\begin{align}
\min_{\boldsymbol \theta,\mathbf{u},\sigma_s^2}\frac{1}{2N}\|\mathbf{y}-\mathbf{H}(\boldsymbol \theta)\mathcal{S}_{\sigma_s}(\mathbf{f}(\mathbf{x};\mathbf{u}))\|^2_2+\lambda \sigma_s^2,\label{metric}
\end{align}

where the first term corresponds to the Mean Squared Error (MSE), which gives an estimate of the observed noise variance, and where the second term is a regularization term that depends on the data noise variance. Note that the network training using the pilot symbols available at the receiver, $\mathbf{x}$, leads to the joint estimation of the channel parameters and symbol detection through the obtained vector $\widehat{\mathbf{u}}$.

The regularization parameter $\lambda$ plays an important role to prevent the overfitting problem. Indeed, for small number of pilot samples $\mathbf{x}$, only minimizing the estimated observed noise variance tends to spread the additive noise into the vector $\mathbf{u}$. This noise contribution could have a negative effect on the Symbol Error Rate (SER), which is our main objective. Increasing $\lambda>0$ limits this phenomenon by encouraging a small data noise variance.

§ SIMULATION RESULTS

This section illustrates the advantages of our approach using a SISO communications system composed of 5 layers: a transmitter IQ imbalance, a FIR channel with $\mathbf{h} =[0.9+0.1j, 0.1+0.1j, 0.01+0.05j,0.02-0.003j, 0.004+0.012j]$, a receiver IQ imbalance, a time-varying Wiener phase noise with instantaneous phase $\varphi[n]=\sum_{k=0}^{n}\psi[n]$ where $\psi[n]\sim \mathcal{N}(0,\sigma_p^2)$, and a Gaussian additive noise $b[n]\sim \mathcal{N}(0,\frac{\sigma_b^2}{2})$. The simulation parameters are provided in Table <ref>. The transmitted signal $\mathbf{s}$ is composed of $N=200$ symbols generated from a 16QAM constellation.

To jointly estimate $\boldsymbol \theta$ and detect $\mathbf{u}$, we consider a multi-layer parametric network composed of 6 layers: a data-augmented layer with $N_u=180$ complex-valued data symbols, a soft projector with initial parameter $\sigma_s^2=1$, a transmitter IQ imbalance layer initialized with $\boldsymbol \theta=[1,0,0,0]$, a FIR layer initialized with $\mathbf{h} =[1,0, 0, 0, 0]$, a receiver IQ imbalance layer initialized with $\boldsymbol \theta=[1,0,0,0]$, and a Piecewise-Constant Phase Layer initialized with $N/N_s=20$ zeros. Therefore, the network is composed of $1+2\times 180+4+2\times 5+4+20=399$ real-valued trainable parameters. The proposed network has been implemented using PyTorch[The source code is available at <https://github.com/vincentchoqueuse/ParamNET>]. The input vector $\mathbf{x}$ is composed of $N_x = 20$ complex-valued symbols transmitted periodically (one pilot symbol followed by 9 unknown data symbols). The loss function is given by (<ref>) and the regularization parameter is set to $\lambda=0.001$. In each experiment, the training stage uses the ADAM optimizer with a learning rate of $10^{-3}$ <cit.>. The performance of the network is evaluated by computing the SER for the complex-valued estimated symbols $\widehat{u}_c[n]=\widehat{u}[n]+j\widehat{u}[n+N_u]$ after projection on the alphabet set $\mathcal{M}$, where $\widehat{u}[n]$ corresponds to the $n^{th}$ trainable parameter of the first layer.

Layer     Type     Parameters

1     IQ Tx     $(\mu,\nu)=(0.9-0.4j,0.4+0.1j)$

2     FIR Channel     $\mathbf{h}$

3     IQ Rx     $(\mu,\nu)=(1.8+0.13j,0.1+0.2j)$

4     Phase Noise     $\sigma_p^2=0.000125$

5     Noise     $\sigma_b^2$

Communication Chain Parameters.

Method     Activation Function     $N_{it}$

Simple     $\times$     $\times$

PG$\_$500     $\times$     500 iterations

PG$\_$1000     $\times$     1000 iterations

PG$\_$2000     $\times$     2000 iterations

Proposed    soft thresholding     $\times$

Considered Techniques

§.§ Single Trial Analysis

Fig. <ref> shows the evolution of $\sigma_s^2$, the MSE and the SER during the training stage. The performance of the proposed technique is compared with the performance obtained with the same parametric network without activation layer. While the MSE is clearly lower for the simple parametric network, the SER shows that this network strongly overfits the data due to the small number of pilots. After 20000 iterations, the SER is equal to $0.0361$ for the parametric network, while the proposed network leads to a SER equal to $0.0027$.

ymin = 0.1,
xlabel=Number of iterations,
xticklabel style=
/pgf/number format/fixed,
/pgf/number format/precision=4
scaled x ticks=false,
[blue,thick] table[x index=0, y index=1,col sep=comma] ./data/one_shot/one_shot.csv;

Evolution of the soft thresholding parameter $\sigma_s^2$

xlabel=Number of iterations,
xticklabel style=
/pgf/number format/fixed,
/pgf/number format/precision=4
xtick=5000,10000,15000, 20000,
scaled x ticks=false,
ylabel=MSE (training)]
[blue,thick] table[x index=0, y index=2,col sep=comma] ./data/one_shot/one_shot.csv;
[red,thick] table[x index=0, y index=3,col sep=comma] ./data/one_shot/one_shot.csv;
No activation, soft thresholding

Evolution of the MSE (training performance)

xlabel=Number of iterations,
xticklabel style=
/pgf/number format/fixed,
/pgf/number format/precision=4
xtick=5000,10000,15000, 20000,
scaled x ticks=false,

ylabel=SER (testing)]
[blue,thick] table[x index=0, y index=4,col sep=comma] ./data/one_shot/one_shot.csv;
[red,thick] table[x index=0, y index=5,col sep=comma] ./data/one_shot/one_shot.csv;
No activation, soft thresholding

Evolution of the Symbol Error Rate (SER) (testing performance)

Evolution of $\sigma_s^2$ and SER performance (SNR=$20$dB).

xlabel=SNR (dB),
legend pos=south west,
[blue,thick] table[x index=0, y index=1,col sep=comma] ./data/monte_carlo/ser.csv;
[thick,brown] table[x index=0, y index=2,col sep=comma] ./data/monte_carlo/ser.csv;
[thick,black] table[x index=0, y index=3,col sep=comma] ./data/monte_carlo/ser.csv;
[thick,magenta,dashed] table[x index=0, y index=4,col sep=comma] ./data/monte_carlo/ser.csv;
[thick,red] table[x index=0, y index=5,col sep=comma] ./data/monte_carlo/ser.csv;
Simple, PG$\_$500 ,PG$\_$1000,PG$\_$2000,Proposed

Testing performance versus SNR.

§.§ Monte Carlo Simulations

Fig. <ref> presents the average SER versus SNR using Monte Carlo simulations. In each Monte Carlo trial, we generate a new realization of the Wiener phase and additive noises, and the optimization is stopped after $N_{it}=10000$ iterations. The performance of the proposed technique is compared with the performances of the techniques reported in Table <ref>.
The Projected Gradient (PG) techniques are obtained by periodically projecting the trainable parameters $\mathbf{u}$ into the set $\mathcal{M}$ after $N_{it}$ iterations <cit.>. As the projection rate has a significant impact on the SER performance, we have considered several PG detectors with different values of $N_{it}$. We observe that the simple technique requires a large SNR to perform well. Furthermore, we also note that the performance of the PG techniques critically depends on the projection rate and SNR. In particular, projecting frequently $\mathbf{u}$ into the set $\mathcal{M}$ leads to poor performance for medium or high SNR. Regarding the proposed method, it leads to the best performance at medium and high SNR.

§ CONCLUSION

This paper proposes a new network architecture called ParamNet for joint estimation and data detection in communications systems. ParamNet is composed of a data augmented layer and a soft thresholding layer that encourages the data constellation. This network can be trained using a small number of pilot symbols and a custom regularized loss function. Simulation results have shown that this strategy outperforms projected gradient techniques. In particular, the proposed strategy provides a solution to avoid the overfitting problem in the medium and high SNR regions. Future works will investigate the application of the proposed network in more complex communications chains, including MIMO systems.