# A Spectral Approach to Optimal Control of the Fokker-Planck Equation

Dante Kalise, Lucas Moschen, Grigorios Pavliotis, Urbain Vaes

► To cite this version:

Dante Kalise, Lucas Moschen, Grigorios Pavliotis, Urbain Vaes. A Spectral Approach to Optimal Control of the Fokker-Planck Equation. IEEE Control Systems Letters, 2025, 9, pp.504-509. 10.1109/LCSYS.2025.3573604. [hal-05453409](https://HAL Id: hal-05453409)

https://inria.hal.science/hal-05453409v1

Submitted on 12 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A Spectral Approach to Optimal Control of the Fokker-Planck Equation*

Published in IEEE Control Systems Letters (L-CSS) at https://ieeexplore.ieee.org/document/11015582

Dante Kalise<fnref n="1" />, Lucas M. Moschen<fnref n="1" />, Grigorios A. Pavliotis<fnref n="1" /> and Urbain Vaes<fnref n="2" />

**Abstract**— In this paper, we present a spectral optimal control framework for Fokker-Planck equations based on the standard ground state transformation that maps the Fokker-Planck operator to a Schrödinger operator. Our primary objective is to accelerate convergence toward the (unique) steady state. To fulfill this objective, a gradient-based iterative algorithm with Pontryagin's maximum principle and the Barzilai-Borwein update is developed to compute time-dependent controls. Numerical experiments on two-dimensional ill-conditioned normal distributions and doublewell potentials demonstrate that our approach effectively targets slow-decaying modes, thus increasing the spectral gap.

## I. INTRODUCTION

Optimal control of diffusion processes, whose dynamics are modeled by Fokker-Planck equations (FPEs), arises in a wide range of applications, from molecular dynamics [1] and optics [2] to social phenomena [3], [4], [5], [6]. By designing controls that act on the FPE, one can steer the probability distribution toward a desired state, accelerate convergence to equilibrium, or stabilize otherwise slow-mixing dynamics. In this work, we focus on the linear FPE as a fundamental case study, noting that the same framework can in principle be extended to nonlinear or mean-field PDEs via linearization.

In controlling the FPE, first- and second-order optimality conditions have been derived for an optimal control problem [7]. Meanwhile, infinite-horizon control strategies for overdamped Langevin dynamics in bounded domains have been proposed to stabilize the FPE and enhance convergence to the stationary measure (Gibbs) [8], [9], both using finite-difference schemes. In particular, the ground-state change of variables we present in **Section II** is the same as [9]. This is especially relevant to sampling, where optimal Langevin

samplers improve mixing by modifying the drift term of the associated SDE [10]. A discussion of the FP control framework, including its formulation for various FPEs and its connection to dynamic programming, can be found in [11].

Unlike finite-volume or finite-difference discretizations, spectral methods rely on eigenfunction expansions to approximate solutions accurately while maintaining low computational cost. This approach is especially advantageous in unbounded domains, where conventional grid-based methods struggle. In addition, spectral formulations naturally connect to Schrödinger-type transformations of the Fokker-Planck operator, enabling direct use of techniques from quantum mechanics and numerical linear algebra [12]. This link was successfully exploited in multiscale stochastic differential equations, demonstrating the versatility of spectral approaches [13].

To our knowledge, no prior work has combined, in an unbounded-domain setting, (i) the ground-state (Schrödinger) transform on $\mathbb{R}^d$, (ii) spectral discretization of the FPE, (iii) open-loop Pontryagin optimal control solved with Barzilai-Borwein (BB) solver, and (iv) systematic shape-function design to control the slowest eigenmodes. The main contributions of our paper are:

* We extend the Hermite-based approach of [14] from quadratic to smooth confining potentials, relying on the spectrum of the associated Schrödinger operator and time-dependent controls.

* We derive the state-adjoint optimality system for a finite-horizon cost and propose a gradient-based iterative algorithm with a BB update to compute the optimal control.

* We extend the bounded-domain shape-function design of [8] to $\mathbb{R}^d$ to target the slow-decaying spectral modes, and we demonstrate via numerical experiments that our method accelerates convergence toward the steady state.

The remainder of this paper is structured as follows: **Section II** provides a theoretical background on the Fokker-Planck equation and the Schrödinger operator. **Section III** formulates the optimal control problem. **Section IV** presents our numerical approach, detailing the spectral discretization. Finally, **Section V** provides numerical experiments.

**Notation.** We denote by $L^2(\mathbb{R}^d)$ the space of square-integrable functions with respect to the Lebesgue measure. For a given probability density $\mu$ on $\mathbb{R}^d$, we say that $\phi$ belongs to the weighted space $L^2(\mathbb{R}^d, \mu)$ if $\phi\sqrt{\mu} \in L^2(\mathbb{R}^d)$. For a given Hilbert space $(H, \langle \cdot, \cdot \rangle)$, we denote $\| \cdot \|_H$ the norm in $H$ defined through the inner product of $H$.

<footnote>*DK is partially supported by the EPSRC Standard Grant EP/T024429/1. LM is supported by the CNRS-ICL PhD Studentship program. GP is partially supported by an ERC-EPSRC Frontier Research Guarantee through Grant No. EP/X038645, ERC Advanced Grant No. 247031 and a Leverhulme Trust Senior Research Fellowship, SRF\R1\241055. UV is partially supported by the European Research Council (ERC) under the EU Horizon 2020 program (grant agreement No 810367) and by the Agence Nationale de la Recherche under grants ANR-21-CE40-0006 (SINEQ) and ANR-23-CE40-0027 (IPSO).</footnote>

<footnote n="1">¹DK, LM, and GP are with the Dept. of Mathematics, Imperial College, London, UK {d.kalise-balza,lmm122,pav}@ic.ac.uk</footnote>

<footnote n="2">²UV is with Inria, Paris, France urbain.vaes@inria.fr</footnote>

<footnote>Published in IEEE Control Systems Letters (L-CSS): https://ieeexplore.ieee.org/document/11015582</footnote>

<footnote>© 2025 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.</footnote>

<footnote>For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising.</footnote>
## II. FOKKER-PLANCK EQUATION AND PRELIMINARIES

Consider the diffusion process $\{X_t\}_{t \ge 0}$ in $\mathbb{R}^d$ that satisfies the stochastic differential equation

$$dX_t = -\nabla V(X_t) dt + \sqrt{2\sigma} dW_t, \quad (1)$$

where $V: \mathbb{R}^d \to \mathbb{R}$ is a potential function, $\sigma > 0$ is the diffusion coefficient, and $\{W_t\}_{t \ge 0}$ denotes a standard $d$-dimensional Brownian motion. The evolution of the probability density function (PDF) $\rho(x, t)$ of $X_t$ is governed by the forward Fokker-Planck equation

$$\partial_t \rho = \mathcal{L}^* \rho, \quad (2)$$

with the initial condition $\rho(x, 0) = \rho_0(x)$, which represents the Lebesgue density function of $X_0$, satisfying $\rho_0 \ge 0$, $\int_{\mathbb{R}^d} \rho_0(x) dx = 1$, and $\rho_0 \in C^\infty(\mathbb{R}^d)$. The operator $\mathcal{L}^*: L^2(\mathbb{R}^d, \rho_\infty^{-1}) \to L^2(\mathbb{R}^d, \rho_\infty^{-1})$ is defined as

$$\mathcal{L}^* \phi := \nabla \cdot (\phi \nabla V) + \sigma \Delta \phi.$$

We assume that $V \in C^\infty(\mathbb{R}^d)$ and that

$$\lim_{|x| \to \infty} V(x) = \infty, \text{ and } e^{-V(\cdot)/\sigma} \in L^1(\mathbb{R}^d)$$

for all $\sigma > 0$. Under these assumptions, the PDE is well-posed (see [14, Thm. 2], [12, Thm. 4.2], [8, Prop. 2.1]) and admits a unique stationary density

$$\rho_\infty(x) := Z^{-1} e^{-V(x)/\sigma},$$

$$\text{where } Z := \int_{\mathbb{R}^d} e^{-V(x)/\sigma} dx \text{ [12, Prop. 4.6].}$$

Next, we define

$$W(x) := \frac{1}{4\sigma} |\nabla V(x)|^2 - \frac{1}{2} \Delta V(x) \quad (3)$$

and assume that $\lim_{|x| \to \infty} W(x) = \infty$. Then $\rho_\infty$ satisfies a Poincaré inequality with a constant $\lambda > 0$ [15, Cor. 1.6]; in particular, if $\rho_0 \in L^2(\mathbb{R}^d, \rho_\infty^{-1})$,

$$\|\rho(\cdot, t) - \rho_\infty\|_{L^2(\mathbb{R}^d, \rho_\infty^{-1})} \le e^{-\lambda \sigma t} \|\rho_0 - \rho_\infty\|_{L^2(\mathbb{R}^d, \rho_\infty^{-1})}.$$

Under these hypotheses, the backward generator $\mathcal{L}$ associated with (1) is self-adjoint in $L^2(\mathbb{R}^d, \rho_\infty)$ and has compact resolvent [12, Sec 4.4].

We perform the (ground state) unitary transformation to map the Fokker-Planck operator $\mathcal{L}^*$ to a Schrödinger operator defined in $L^2(\mathbb{R}^d)$. Specifically, define the unitary operator

$$\mathcal{U}: L^2(\mathbb{R}^d, \rho_\infty^{-1}) \to L^2(\mathbb{R}^d) \text{ by } (\mathcal{U}\phi)(x) := \frac{\phi(x)}{\sqrt{\rho_\infty(x)}}.$$

Verifying that $\mathcal{U}$ is an isometric isomorphism with a well-defined inverse is straightforward.

Based on the unitary operator, we define a new operator $\mathcal{H}$ on $L^2(\mathbb{R}^d)$ by

$$\mathcal{H} := -\mathcal{U} \mathcal{L}^* \mathcal{U}^{-1},$$

which has the form $\mathcal{H} = -\sigma\Delta + W(x)$. The operator $\mathcal{H}$ is self-adjoint in $L^2(\mathbb{R}^d)$ and possesses a purely discrete, unbounded spectrum with eigenvalues of finite multiplicity [16, Thm. XIII.67]. Moreover, standard elliptic regularity theory guarantees that its eigenfunctions are infinitely differentiable.

Under the transformation $\psi := \mathcal{U}\rho$, the Fokker-Planck PDE transforms to the equation

$$\partial_t \psi = -\mathcal{H}\psi, \quad (4)$$

with initial condition $\psi(\cdot, 0) = \psi_0 := \mathcal{U}\rho_0$. Equation (4) resembles a Schrödinger equation with imaginary time. By transforming the Fokker-Planck operator into a Schrödinger operator, we recast the problem in the standard flat $L^2(\mathbb{R}^d)$ space, where the resulting operator is self-adjoint and its spectral properties can be analyzed using well-established results from the theory of Schrödinger operators.

## III. OPTIMAL CONTROL PROBLEM FORMULATION

We consider that $X_t$ is subject to an additional conservative force:

$$dX_t = -\nabla V(X_t) dt - \sum_{i=1}^{m} u_i(t) \nabla \alpha_i(X_t) dt + \sqrt{2\sigma} dW_t.$$

Here, $u_i: [0, T] \to \mathbb{R}$ are time-dependent controls in $L^2([0, T])$, and $\alpha_i: \mathbb{R}^d \to \mathbb{R}$ are shape control functions in $C^\infty(\mathbb{R}^d)$. In practice, $-\nabla \alpha_i$ can be generated by optical near-field traps [17] or by virtual guidance in crowd and robotic applications. The associated FPE is given by

$$\partial_t \rho = \mathcal{L}^* \rho + \sum_{i=1}^{m} u_i(t) A_i^* \rho, \quad (5)$$

$$\text{where, for } i = 1, \dots, m, A_i^* \phi := \nabla \cdot (\phi \nabla \alpha_i(x)).$$

The objective is to determine a pair $(\rho, \mathbf{u})$ that minimizes the cost functional

$$J(\rho, \mathbf{u}) := \frac{1}{2} \|\rho(\cdot, T) - \rho^\dagger\|_{L^2(\mathbb{R}^d, \rho_\infty^{-1})}^2 + \frac{\nu}{2} \int_0^T |\mathbf{u}(t)|^2 dt + \frac{\kappa}{2} \int_0^T \|\rho(\cdot, t) - \hat{\rho}\|_{L^2(\mathbb{R}^d, \rho_\infty^{-1})}^2 dt \quad (6)$$

subject to the constraint that $\rho$ satisfies equation (5), with $\mathbf{u}(t) \in \mathbb{R}^m$, $T > 0$, $\nu > 0$, $\kappa > 0$, and $\rho^\dagger \in L^2(\mathbb{R}^d, \rho_\infty^{-1})$ is a desired final distribution. We also include a probability distribution $\hat{\rho}$. We focus on the case where $\rho^\dagger = \hat{\rho} = \rho_\infty$.

Applying the unitary transformation $\mathcal{U}$ and defining $\mathcal{N}_i = \mathcal{U} A_i^* \mathcal{U}^{-1}$, we transform the system into the controlled (imaginary-time) Schrödinger equation. The problem then reduces to finding a pair $(\psi, \mathbf{u})$ that minimizes

$$J(\psi, \mathbf{u}) := \frac{1}{2} \|\psi(\cdot, T) - \psi^\dagger\|_{L^2(\mathbb{R}^d)}^2 + \frac{\nu}{2} \int_0^T |\mathbf{u}(t)|^2 dt + \frac{\kappa}{2} \int_0^T \|\psi(\cdot, t) - \hat{\psi}\|_{L^2(\mathbb{R}^d)}^2 dt, \quad (7)$$

where $\psi^\dagger := \mathcal{U}\rho^\dagger$ and $\hat{\psi} := \mathcal{U}\hat{\rho}$, subject to the controlled equation

$$\partial_t \psi = -\mathcal{H}\psi + \sum_{i=1}^{m} u_i(t) \mathcal{N}_i \psi, \quad (8)$$

with each $\mathcal{N}_i$ modulating $\psi$ via $u_i(t)$ (see Proposition 1 for its explicit form).

**Proposition 1:** Define $b_i := \nabla \alpha_i \cdot \nabla \log \sqrt{\rho_\infty}$. Then, for all $\phi \in L^2(\mathbb{R}^d)$,

$$\mathcal{N}_i \phi = \nabla \cdot (\phi \nabla \alpha_i) + b_i \phi.$$
Moreover, the formal adjoint of $\mathcal{N}_i$ in $L^2(\mathbb{R}^d)$ is

$$ \mathcal{N}_i^* \phi = b_i \phi - \nabla \alpha_i \cdot \nabla \phi. $$

*Proof:* See the Appendix (page 6).

<figure><img src="image_2.png" /></figure>

Standard results in infinite-dimensional optimal control (see, e.g., [18, Ch. 3]) guarantee the existence of an optimal control given the coercivity and weak lower semicontinuity of the cost functional, together with the continuity of the control-to-state mapping. In our setting, the quadratic structure of the cost and the compactness of the resolvent of $\mathcal{H}$ ensure that a minimizing pair $(\psi^*, \mathbf{u}^*)$ exists in the admissible control space $L^2(0, T; \mathbb{R}^m)$. A similar result can be found in [19, Thm. 1]. In this reference, the Schrödinger equation from quantum mechanics (i.e. in “real-time”) is considered. However, the proof is based on the spectral analysis of the elliptic operator and it can be applied directly to our problem.

To solve problem (7)–(8), we consider the first-order optimality conditions. These conditions are derived by introducing a function $\varphi \in L^2(\mathbb{R}^d)$ as the adjoint variable and forming the augmented Lagrangian

$$ \mathcal{L}(\psi, \mathbf{u}, \varphi) := J(\psi, \mathbf{u}) + \int_0^T \langle \varphi(\cdot, t), \partial_t \psi(\cdot, t) \rangle_{L^2(\mathbb{R}^d)} dt \\ + \int_0^T \langle \varphi(\cdot, t), \mathcal{H}\psi(\cdot, t) - \sum_{i=1}^m u_i(t) \mathcal{N}_i \psi(\cdot, t) \rangle_{L^2(\mathbb{R}^d)} dt. $$

The first-order optimality conditions are then obtained by taking variations with respect to $\psi$, $\mathbf{u}$, and $\varphi$, and setting the corresponding gradients to zero.

The adjoint variable $\varphi$ satisfies

$$ \left\{ \begin{array}{l} \partial_t \varphi = \mathcal{H}\varphi - \sum_{i=1}^m u_i(t) \mathcal{N}_i^* \varphi + \kappa(\psi - \hat{\psi}), \\ \varphi(x, T) = \psi^\dagger(x) - \psi(x, T), \end{array} \right. \qquad (9) $$

with the optimal control given by

$$ \nabla_{u_i} \mathcal{L} := \nu u_i - \langle \mathcal{N}_i^* \varphi, \psi \rangle_{L^2(\mathbb{R}^d)} = 0, \qquad (10) $$

$$ u_i(t) = \frac{1}{\nu} \langle \varphi, \mathcal{N}_i \psi \rangle_{L^2(\mathbb{R}^d)} = \frac{1}{\nu} \langle \mathcal{N}_i^* \varphi, \psi \rangle_{L^2(\mathbb{R}^d)}. \qquad (11) $$

In equation (9), we used that $\mathcal{H}$ is a self-adjoint operator. Therefore, an optimal pair $(\psi^*, \mathbf{u}^*)$ satisfies equations (8), (9), and (11).

## IV. NUMERICAL METHOD: SPECTRAL DISCRETIZATION

Let $\{\lambda_k, e_k\}_{k \ge 0}$ denote the eigenvalues and eigenfunctions of the operator $\mathcal{H}$. Without loss of generality, we reorder the eigenvalues as a nondecreasing sequence so that $\lambda_0 = 0$ and $\lambda_1 > 0$. The eigenfunctions $\{e_k\}_{k \ge 0}$ form a $L^2(\mathbb{R}^d)$-orthonormal basis with $e_0 = \sqrt{\rho_\infty}$ since $\rho_\infty$ is the steady state of equation (2).

The solution $\psi(x, t)$ of equation (8) is expanded as

$$ \psi(x, t) = \sum_{k=0}^{\infty} a_k(t) e_k(x). $$

Substituting this expansion into the equation (8) and projecting onto each $e_i(x)$ yields the infinite system of ordinary differential equations (ODEs)

$$ \dot{a}_i(t) = -\lambda_i a_i(t) + \sum_{j=1}^{m} u_j(t) \sum_{k=0}^{\infty} a_k(t) \langle \mathcal{N}_j e_k, e_i \rangle_{L^2(\mathbb{R}^d)} \quad (12) $$

Since $\lambda_k \to \infty$ as $k \to \infty$, higher modes decay rapidly, which motivates the truncation of the expansion to a finite-dimensional system.

**Remark.** By construction $a_0(t) = \langle \psi(\cdot, t), e_0 \rangle_{L^2(\mathbb{R}^d)} \equiv 1$, implying that $\dot{a}_0(t) = 0$, which agrees with $\mathcal{N}_i^*[\sqrt{\rho_\infty}] = 0$.

Similarly, expanding $\varphi(x, t) = \sum_{k=0}^\infty p_k(t) e_k(x)$ leads to

$$ \begin{aligned} \dot{p}_i(t) = \lambda_i p_i(t) & - \sum_{j=1}^{m} u_j(t) \sum_{k=0}^{\infty} p_k(t) \langle \mathcal{N}_j^* e_k, e_i \rangle_{L^2(\mathbb{R}^d)} \quad (13) \\ & + \kappa(a_i - \hat{a}_i), \end{aligned} $$

where $\hat{a}_i$ are the spectral coefficients corresponding to $\hat{\psi}$. The optimal control is then given by

$$ u_i(t) = \frac{1}{\nu} \sum_{j,k=0}^{\infty} a_j(t) p_k(t) \langle \mathcal{N}_i^* e_k, e_j \rangle_{L^2(\mathbb{R}^d)}. \qquad (14) $$

Considering a truncation to $N$ eigenmodes, we define the matrices $B^i, A^i \in \mathbb{R}^{N \times N}$ by

$$ B_{jk}^i := \langle b_i e_j, e_k \rangle_{L^2(\mathbb{R}^d)}, \quad A_{jk}^i := \langle \nabla \alpha_i \cdot \nabla e_j, e_k \rangle_{L^2(\mathbb{R}^d)}. $$

The resulting coupled finite-dimensional system is

$$ \left\{ \begin{aligned} \dot{\mathbf{a}}(t) &= -\Lambda \mathbf{a}(t) + \sum_{j=1}^{m} u_j(t) (B^j - A^j) \mathbf{a}(t), \\ \dot{\mathbf{p}}(t) &= \Lambda \mathbf{p}(t) - \sum_{j=1}^{m} u_j(t) (B^j - (A^j)^T) \mathbf{p}(t) \\ &\quad + \kappa(\mathbf{a}(t) - \hat{\mathbf{a}}), \end{aligned} \right. \qquad (15) $$

with initial and terminal conditions

$$ a_i(0) = \langle \psi_0, e_i \rangle_{L^2(\mathbb{R}^d)}, \quad p_i(T) = \langle \psi^\dagger, e_i \rangle_{L^2(\mathbb{R}^d)} - a_i(T), $$

where $\Lambda$ is a diagonal matrix such that $\Lambda_{ii} = \lambda_{i-1}$ for $i \in \{1, \dots, N\}$ and $B^j$ matrices are symmetric.

**Remark.** When $\psi^\dagger = \hat{\psi} = \sqrt{\rho_\infty}$, i.e., the target state is the steady state, convergence is achieved provided that $\sum_{i=1}^N a_i(t)^2 \to 0$, i.e. when all modes apart from the steady-state mode decay to zero. Since modes corresponding to small eigenvalues $\lambda_i$ decay more slowly, it is important to design the control to target and accelerate the decay of them.

Under the regularity assumptions on $W$, classical spectral convergence theory ensures that the projection error decays rapidly with the number of modes. In particular, the truncated expansion converges to $\psi$ in the $L^2$ norm at a rate that depends on the smoothness of $\psi$ (see, e.g., [14]).

### A. The Iterative Method

To compute the optimal control numerically, we employ a gradient descent method applied to the coupled system (15). The space-discretization of the gradient (10) is given by

$$ \nabla_{u_j} J := \nu u_j - \langle \mathbf{p}(t), (B^j - A^j)\mathbf{a}(t) \rangle_{\mathbb{R}^N}. \qquad (16) $$

Using this gradient, we perform the following iterative procedure, outlined in Algorithm 1.

In our numerical implementation, the eigenfunctions $\{e_k\}$ and eigenvalues $\{\lambda_k\}$ of the operator $\mathcal{H}$ are computed using a finite element approach implemented in Wolfram Mathematica via the built-in function NDEigensystem [20] under Dirichlet boundary conditions on a sufficiently large domain, ensuring that the potential $W$ is large. We emphasize that
the eigenpair computation is treated as a black box: while we use finite elements in 2D for its robustness, any accurate eigensolver, such as global spectral discretizations, low-rank tensor methods, or physics-informed neural networks [21], can be used without affecting our convergence results.

Once the eigenfunctions are computed, we use a trapezoidal rule to obtain the matrices $B^i$ and $A^i$ for $i \in \{1, \dots, m\}$. In a full-grid (tensor-product) spectral discretization retaining $N$ modes per coordinate, the total number of basis functions grows as $N^d$ and each coupling matrix $A^j, B^j$ has size $N^d \times N^d$, so memory and computational costs scale like $O(N^{2d})$.

In Algorithm 1, the forward and adjoint systems are integrated using a Runge-Kutta method 5(4) with relative tolerance $1 \times 10^{-7}$ and absolute tolerance $1 \times 10^{-9}$. The BB step size is computed following the formulation as in [22]. For strictly convex quadratic objectives, BB enjoys R-linear convergence, i.e. exponential decay of the error toward the global minimizer, and even Q-linear convergence, meaning each iteration shrinks the distance to the minimizer by a fixed factor less than one, whenever the Hessian's condition number satisfies $\kappa < 2$ [22], [23]. For nonconvex objectives with Lipschitz gradients and bounded level sets, BB still drives the gradient to zero [24]. In our spectral-PMP problem, the reduced cost $J$ is smooth, coercive, strongly convex in $u$, and has bounded sublevel sets. Hence, $\|\nabla J(u^k)\| \to 0$ even though the iterates need not decrease monotonically. Other approaches, such as line searches with Wolfe conditions or the nonlinear conjugate-gradient with the Armijo condition, can be employed, but we chose the BB update for its simplicity and efficiency.

**Algorithm 1** Spectral control solver with reduced gradient and Barzilai-Borwein update

**Require:** $\text{tol} > 0$, $k_{\text{max}}$, initial control function $u^0(t)$ (and $u^{-1}(t)$ for BB), final time $T$, and a time grid $t_{\text{eval}}$.

1: **Initialize:** $k \leftarrow 0$.

2: **while** $\|\nabla J(u^k)\| > \text{tol}$ and $k < k_{\text{max}}$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;**(1)** Compute $\mathbf{a}(t)$ at $t_{\text{eval}}$ by solving (12) with control $\mathbf{u}^k$ and initial condition $\mathbf{a}(0)$.

4: &nbsp;&nbsp;&nbsp;&nbsp;**(2)** Compute $\mathbf{p}(t)$ at $t_{\text{eval}}$ by solving (13) with controls $\mathbf{u}^k$, state $\mathbf{a}(t)$, and terminal condition $\mathbf{p}(T) = \mathbf{a}^\dagger - \mathbf{a}(T)$.

5: &nbsp;&nbsp;&nbsp;&nbsp;**(3)** Evaluate $\nabla J_{u_j}(u_j^k)$ from (16) for each $j$.

6: &nbsp;&nbsp;&nbsp;&nbsp;**(4)** Compute the Barzilai-Borwein step size

$$ \gamma_k = \frac{\sum_{j=1}^{m} \langle u_j^k - u_{j-1}^k, \nabla J_{u_j}(u_j^k) - \nabla J_{u_j}(u_{j-1}^k) \rangle}{\sum_{j=1}^{m} |\nabla J_{u_j}(u_j^k) - \nabla J_{u_j}(u_{j-1}^k)|^2} $$

7: &nbsp;&nbsp;&nbsp;&nbsp;**(5)** For each $j \in \{1, \dots, m\}$, update the control:

$$ u_j^{k+1}(t) = u_j^k(t) - \gamma_k \nabla J_{u_j}(u_j^k(t)). $$

8: &nbsp;&nbsp;&nbsp;&nbsp;**(6)** Update the control interpolants on $t_{\text{eval}}$.

9: &nbsp;&nbsp;&nbsp;&nbsp;**(7)** Set $k \leftarrow k + 1$.

10: **end while**

11: **return** Optimal control $\mathbf{u}^k(t)$, state coefficients $\mathbf{a}(t)$, and adjoint coefficients $\mathbf{p}(t)$.

## B. Choice of the Shape Function

In our optimal control framework, we only optimize over the time-dependent controls $u_i(t)$ while keeping the shape functions $\alpha_i(x)$ fixed. These shape functions modulate the spatial distribution of the control input and influence the spectral properties of the operator $\mathcal{H} - u_i N_i$. In particular, a careful choice of $\alpha_i$ can increase the spectral gap of the controlled operator, thereby accelerating the convergence of the system toward its steady state.

For the purpose of choosing the functions $\alpha_i$, we introduce the deviation variable $y(x, t) = \psi(x, t) - \sqrt{\rho_\infty(x)}$ for every $(x, t) \in \mathbb{R}^d \times \mathbb{R}_{\ge 0}$. The deviation function $y(x, t)$ satisfies

$$ y_t = -\mathcal{H}y + \sum_{j=1}^{m} u_j N_j y + \sum_{j=1}^{m} u_j B_j, \quad (17) $$

where $B_j := N_j[\sqrt{\rho_\infty}]$ represents the effect of the control on the steady state.

Motivated by the Lyapunov-based feedback strategy introduced in [8], we select the shape function $\alpha_j \in L^2(\mathbb{R}^d)$ to align the control with the slowest decaying modes. Specifically, we choose $\alpha_j$ such that

$$ B_j = e_j, \quad \text{i.e.} \quad \nabla \cdot (\rho_\infty \nabla \alpha_j) = \sqrt{\rho_\infty} e_j, $$

for $j \in \{1, \dots, m\}$ the $m$ slowest modes of the operator $\mathcal{H}$. This condition intuitively ensures that the control acts in the direction of the slowest decaying modes, which are typically the bottleneck for convergence. The above PDE has a unique solution in $L^2(\mathbb{R}^d)$. Moreover, by the smoothness of $\rho_\infty$ and $e_j$, the solution $\alpha_j$ is smooth as well [25].

Although we do not claim that this choice is optimal in any formal sense, it provides a practical means of enhancing the influence of the control on the dynamics. In practice, we approximate the solution of this PDE by projecting on a finite-dimensional basis, that is, $\alpha_j(x) \approx \sum_{k=1}^N c_{j,k} e_k(x)$, and then solve the resulting linear system for $c_{j,k}$.

## C. Initialization of the Control Functions

Due to the inherent nonconvexity of the optimal control problem – both in its infinite-dimensional formulation and its finite-dimensional discretization – gradient-based methods are not guaranteed to find the global minimum. Such methods may become trapped in suboptimal local minima and exhibit slow convergence. In addition, initializing the control functions as constant zero functions (i.e., $u_i \equiv 0$ for all $i$) may yield a terminal state that is close to the desired one for some parameter choices, but not necessarily accelerating convergence. This behavior can limit the efficiency of the algorithm and not steer the iterative process toward the desired global minimum. Therefore, a well-informed initial guess for $u_i$ shortens the optimization time and reduces the risk of converging to a suboptimal local minimum.

Our strategy for control initialization is based on the linearization of the system around the steady state. Specifically, we consider the linearization of equation (17)

$$ y_t = -\mathcal{H}y + \sum_{j=1}^{m} u_j(t) B_j $$
and formulate an infinite-horizon cost functional

$$J_{\infty}(y, \mathbf{u}) := \frac{1}{2} \int_{0}^{\infty} \left[ \kappa \| y(\cdot, t) \|_{L^2(\mathbb{R}^d)}^2 + \nu |\mathbf{u}(t)|^2 \right] dt,$$

which admits an optimal solution of the form

$$u_j^*(t) = -\frac{1}{\nu} \langle \Pi y(\cdot, t), \mathcal{B}_j^* \rangle_{L^2(\mathbb{R}^d)},$$

where $\Pi$ is the solution of the associated Riccati equation

$$-\mathcal{H}\Pi - \Pi\mathcal{H} - \frac{1}{\nu}\Pi\left(\sum_{j=1}^{m} \mathcal{B}_j^* \mathcal{B}_j\right)\Pi + \kappa \mathcal{I} = 0.$$

Here, $\mathcal{I}$ is the identity operator on $L^2(\mathbb{R}^d)$. This Riccati equation admits a unique solution, providing a well-defined optimal feedback law.

We compute the feedback control $u_j^*(t)$ using spectral discretization as in the previous sections. Since we built matrices $\Lambda$, $B^j$, and $A^j$ for the algorithm, we can use them to solve the Ricatti equation, which is reasonably fast. This computed feedback is then used as an initial guess for the iterative gradient-based optimization algorithm from Section IV-A. Our numerical experiments indicate that this initialization accelerates convergence, reduces numerical approximation errors, and helps avoid convergence to suboptimal solutions.

## V. NUMERICAL EXPERIMENTS AND RESULTS

In this section, we illustrate the performance of the proposed spectral control solver (Algorithm 1) on two benchmark cases. For simplicity, we fix $\sigma = 1$ in all experiments. The first example is the quadratic potential

$$V_1(x, y) = \frac{1}{2}(ax^2 + by^2).$$

The invariant measure for this potential is a Gaussian measure with an ill-conditioned (when $b \ll a$) covariance matrix. This potential allows the analytical computation of the eigenvalues and eigenfunctions of the operator $\mathcal{H}$, using the Metafune-Palara-Priola theorem [10]. In fact, for nonnegative integers $m, n$, the eigenvalues are given by $\lambda_{m,n} = an+bm$. Since the spectral gap is $\min(a, b)$, choosing a small value for either $a$ or $b$ results in slow convergence toward the steady state. We refer to such cases as ill-conditioned. The second example is the separable double-well potential

$$V_2(x, y) = (x^2 - c_x)^2 + (y^2 - c_y)^2.$$

This double-well potential is metastable for small values of $\sigma$ (low temperature), meaning probability mass remains trapped near a local minimum for long times before switching.

In all experiments, we set the final time $T = 5$, control regularization parameter $\nu = 10^{-4}$, and running cost regularization parameter $\kappa = 5$. The solution is approximated using $N = 50$ eigenmodes, and the initial condition is chosen as a truncated Gaussian distribution centered at $(-0.2, 0.5)$. For reference, in the quadratic ill-conditioned Gaussian test, the 50th eigenvalue of $\mathcal{H}$ is $\lambda_{50} = 7$, while in the double-well example it is $\lambda_{50} \approx 48$. We ran 500 iterations for the gradient method, by which time the norm of the gradient falls to $10^{-4} - 10^{-3}$. All experiments were performed on a Mac

with an Apple M3 processor (8 cores) and 16 GB of memory. As a baseline in our plots, we construct an infinite-horizon linear-quadratic-regulator LQR, Riccati-based, feedback law via the method introduced in Section IV-C.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Time evolution of the error norm for the quadratic potential with $a = 1$; line styles distinguish values of $b$ (solid for $b = 0.1$, dashed for $b = 0.05$), while colors indicate the control strategy (black: uncontrolled, blue: LQR feedback control, red: optimal control).</figcaption></figure>

Figure 1 illustrates the evolution of the error norm $|a(t) - a^\dagger|_2$, over time for the quadratic potential with $a = 1$ and two choices of $b$ ($b = 0.05$ and $b = 0.1$). The plot compares the optimal control (red) with the uncontrolled case (black) and the LQR feedback control (blue), demonstrating that the optimal control strategy accelerates convergence even in the challenging case of $b = 0.05$. We notice a reduction of two orders of magnitude for both experiments. Four control functions were set to be optimized, which are depicted in Figure 2. The control signals are strongly active during the initial moments of the experiment, after which they decay as the system approaches the steady state, and the first four modes approach zero.

<figure><img src="image_4.png" /><figcaption>Fig. 2. The four optimized control functions for the quadratic potential experiments with $a = 1$ and $b = 0.1$; control is shown only up to time $t = 0.5$ as it remains close to zero for $t > 0.5$.</figcaption></figure>

To assess robustness to imperfect actuator shapes, we perturbed each learned shape function by adding a Gaussian perturbation $r_j$, i.e. $\tilde{\alpha}_j(x) = \alpha_j(x) + \varepsilon r_j(x)$. We then measured the relative increase in final-time error $\|\mathbf{a}(T) - \mathbf{a}^\dagger\|_2$ for $\varepsilon \in \{1, 5, 10\}$ on the quadratic test ($a = 1, b = 0.1$). The resulting percentage increases are reported in Table I. Here, $\varepsilon$ is the $L^1$-norm of the misalignment added to each actuator. Despite a 9.5% increase in final-time error for a moderate misalignment ($\varepsilon = 1$), the method remains largely effective; only large perturbations cause significant performance loss.

Figure 3 presents the time evolution of the same error for the double-well potential. In the uncontrolled cases, the error norm slowly decays because the system becomes temporarily trapped in one of the wells. In contrast, the optimal control
<table><caption>TABLE I
IMPACT OF SHAPE-FUNCTION MISALIGNMENT ON FINAL-TIME ERROR
(QUADRATIC TEST WITH b = 0.1).</caption><thead><tr><th></th><th>ε = 1</th><th>ε = 5</th><th>ε = 10</th></tr></thead><tbody><tr><th>Δ(ε) [%]</th><td>9.5</td><td>17.4</td><td>950</td></tr></tbody></table>

strategy significantly accelerates convergence, as evidenced by a steeper decay of the error, finishing with a final approximation error of three orders of magnitude less. The LQR feedback control was included as a baseline. The two control functions optimized for this example look similar to those of Figure 2. These results confirm that our spectral control approach is robust even for multimodal potentials. Figure 4 compares the evolution of the coefficients under uncontrolled and controlled dynamics. In particular, the second and third components, $a_1$ and $a_2$, are affected by the control, whereas others (e.g., $a_5$) remain essentially unchanged.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Time evolution of the error norm $||a(t) - a^{\dagger}||_2$ for the double-well potential. The plot compares uncontrolled dynamics (dashed), LQR feedback control (dotted), and optimized open-loop control (solid).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Time evolution of the coefficients of the spectral decomposition of the double-well potential.</figcaption></figure>

## VI. CONCLUSIONS

Our numerical experiments indicate that the proposed spectral control solver effectively accelerates convergence to the desired steady state for both quadratic and double-well potentials by targeting the slower-decaying modes associated with the small eigenvalues. Moreover, the same framework can steer an initial distribution $ρ_0$ toward a nonequilibrium target $ρ^\dagger$ (e.g. in multimodal landscapes) by choosing each $α_i(x)$ to align with the dominant eigenfunctions of the linearized error operator $ψ - ψ^\dagger$. However, because we only use time-dependent scalars $u_i(t)$ multiplying fixed spatial profiles $α_i(x)$, exact controllability to $ρ^\dagger$ is not guaranteed. Finally,

as the spatial dimension increases, the required number of spectral modes grows dramatically, and the storage and manipulation of $A^j$ and $B^j$ become infeasible. Promising remedies include low-rank tensor decompositions, model reduction [9], sparse-grid spectral methods, or data-driven modes. The study of these extensions is left for future work.

## APPENDIX

We now prove Proposition 1.

*Proof:* Applying the product rule to $N_i$ yields

$$\nabla \cdot (\sqrt{\rho_{\infty}} \phi \nabla \alpha_i) = \sqrt{\rho_{\infty}} \nabla \cdot (\phi \nabla \alpha_i) + \phi \nabla \sqrt{\rho_{\infty}} \cdot \nabla \alpha_i.$$

Dividing by $\sqrt{\rho_{\infty}}$ and noting that

$$\frac{\nabla \sqrt{\rho_{\infty}}}{\sqrt{\rho_{\infty}}} = \nabla \log \sqrt{\rho_{\infty}},$$

we obtain $N_i \phi = \nabla \cdot (\phi \nabla \alpha_i) + \phi \nabla \alpha_i \cdot \nabla \log \sqrt{\rho_{\infty}}$. This shows the desired expression. The formula for the formal adjoint follows from standard integration by parts. ■

## REFERENCES

[1] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids. Oxford University Press, 2017.

[2] P. H. Jones, O. M. Maragò, and G. Volpe, Optical Tweezers: Principles and Applications. Cambridge: Cambridge University Press, 2015.

[3] D. Helbing, A. Johansson, J. Mathiesen, M. H. Jensen, and A. Hansen, “Analytical approach to continuous and intermittent bottleneck flows,” Phys. Rev. Lett., vol. 97, no. 16, p. 168001, 2006.

[4] C. Castellano, S. Fortunato, and V. Loreto, “Statistical physics of social dynamics,” Rev. Mod. Phys., vol. 81, no. 2, pp. 591–646, 2009.

[5] G. Albi, Y.-P. Choi, M. Fornasier, and D. Kalise, “Mean field control hierarchy,” Appl. Math. Optim., vol. 76, pp. 93–135, 2017.

[6] L. Pareschi and G. Toscani, “Self-similarity and power-like tails in nonconservative kinetic models,” J. Stat. Phys., vol. 124, pp. 747–779, 2006.

[7] M. S. Aronna and F. Troeltzsch, “First and second order optimality conditions for the control of Fokker-Planck equations,” ESAIM Control Optim. Calc. Var., vol. 27, p. 15, 2021.

[8] T. Breiten, K. Kunisch, and L. Pfeiffer, “Control strategies for the Fokker-Planck equation,” ESAIM Control Optim. Calc. Var., vol. 24, no. 2, pp. 741–763, 2018.

[9] ——, “A reduction method for Riccati-based control of the Fokker-Planck equation,” IFAC-PapersOnLine, vol. 50, no. 1, pp. 1631–1636, 2017.

[10] T. Lelievre, F. Nier, and G. A. Pavliotis, “Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion,” J. Stat. Phys., vol. 152, no. 2, pp. 237–274, 2013.

[11] M. Annunziato and A. Borzi, “A Fokker-Planck control framework for stochastic systems,” EMS Surv. Math. Sci., vol. 5, no. 1-2, pp. 65–98, 2018.

[12] G. A. Pavliotis, Stochastic Processes and Applications, ser. Texts in Applied Mathematics. Springer, 2014, vol. 60.

[13] A. Abdulle, G. A. Pavliotis, and U. Vaes, “Spectral methods for multiscale stochastic differential equations,” SIAM/ASA J. Uncertain. Quantif., vol. 5, no. 1, pp. 720–761, 2017.

[14] M. Mohammadi and A. Borzi, “A Hermite spectral method for a Fokker-Planck optimal control problem in an unbounded domain,” Int. J. Uncertainty Quantif., vol. 5, no. 3, 2015.

[15] D. Bakry, F. Barthe, P. Cattiaux, and A. Guillin, “A simple proof of the Poincaré inequality for a large class of probability measures including the log-concave case,” Electron. Commun. Probab., vol. 13, pp. 60–66, 2008.

[16] M. Reed and B. Simon, IV: Analysis of Operators. Elsevier, 1978, vol. 4.

[17] M. A. Zaman, P. Padhy, and L. Hesselink, “Fokker-planck analysis of optical near-field traps,” Scientific Reports, vol. 9, no. 1, p. 9557, 2019.

[18] J. L. Lions, Optimal Control of Systems Governed by Partial Differential Equations. Springer, 1971, vol. 170.
[19] A. P. Peirce, M. A. Dahleh, and H. Rabitz, "Optimal control of quantum-mechanical systems: Existence, numerical approximation, and applications," *Phys. Rev. A*, vol. 37, no. 12, p. 4950, 1988.

[20] W. R. Inc., "Mathematica, Version 14.2," champaign, IL, 2024. [Online]. Available: https://www.wolfram.com/mathematica

[21] H. Jin, M. Mattheakis, and P. Protopapas, "Physics-informed neural networks for quantum eigenvalue problems," in *2022 Int. Joint Conf. Neural Netw. (IJCNN)*, IEEE, 2022, pp. 1–8.

[22] J. Barzilai and J. M. Borwein, "Two-point step size gradient methods," *IMA J. Numer. Anal.*, vol. 8, no. 1, pp. 141–148, 1988.

[23] M. Raydan, "On the Barzilai and Borwein choice of steplength for the gradient method," *IMA J. Numer. Anal.*, vol. 13, no. 3, pp. 321–326, 1993.

[24] B. Azmi and K. Kunisch, "Analysis of the Barzilai-Borwein step-sizes for problems in Hilbert spaces," *J. Optim. Theory Appl.*, vol. 185, no. 3, pp. 819–844, 2020.

[25] D. Gilbarg and N. S. Trudinger, *Elliptic partial differential equations of second order*. Springer US, 1997, pp. 213–267.