finite in the expansion (5.16), we may consider $(\bar{k},\bar{j},\bar{i})\in
H_{2}$ defined in (5.17) such that the term corresponding to this parameter
dominates all the others. Following what we wrote in Case 2 of the proof of
item (i) of Proposition 5.1 above, we may just take $\tau=0$ to discuss
codominance issues. Using (5.16) and (5.27), this term reads as follows, when
$\tau=0$:

$\displaystyle{s_{A,n}}^{\bar{i}}e^{(1-\frac{\bar{k}}{2}){s_{A,n}}}A^{\bar{k}-1}(\bar{k}-\bar{j})\tilde{C}_{\bar{k},\bar{j},\bar{i}}B_{n}^{\bar{k}-1-\bar{j}}$
(5.33) $\displaystyle=$ $\displaystyle
e^{(1-\frac{\bar{k}}{2}-\theta(\bar{k}-1-\bar{j})){s_{A,n}}}{s_{A,n}}^{\bar{i}+\alpha(\bar{k}-1-\bar{j})}A^{\bar{k}-1}(\bar{k}-\bar{j})\tilde{C}_{\bar{k},\bar{j},\bar{i}}\psi_{n}^{\bar{k}-1-\bar{j}}.$

Let us then characterize the set $\bar{E}$ of all $(k,j,i)\in H_{2}$ such that
the corresponding term in (5.16) is of the same size as the term corresponding
to $(\bar{k},\bar{j},\bar{i})$. From (5.33) and (5.30), we see that we need to
have

$1-\frac{\bar{k}}{2}-\theta(\bar{k}-1-\bar{j})=1-\frac{k}{2}-\theta(k-1-j)\mbox{
and }\bar{i}+\alpha(\bar{k}-1-\bar{j})=i+\alpha(k-1-j).$ (5.34)

This way, we keep only the dominant terms in (5.16), and write

$\displaystyle
v_{b,1,0}({s_{A,n}}+\tau)=e^{\frac{\tau}{2}}\left(\sum_{(k,j,i)\in\bar{E}}\right.$
$\displaystyle({s_{A,n}}+\tau)^{i}e^{(1-\frac{k}{2}){s_{A,n}}}A^{k-1}(k-j)\tilde{C}_{k,j,i}B_{n}^{k-1-j}$
$\displaystyle\left.+O({s_{A,n}}^{\gamma}e^{(\frac{1-M}{2}){s_{A,n}}})\right).$

Using (5.33) and (5.34), together with (5.30), we write for all
$\tau\in[0,(2M-2){s_{A,n}}]$ (note that this interval refers to our discussion
in Case 2 of the proof of item (i) of Proposition 5.1):

$\displaystyle v_{b,1,0}({s_{A,n}}+\tau)=$ $\displaystyle
e^{\frac{\tau}{2}}e^{(1-\frac{\bar{k}}{2}-\theta(\bar{k}-1-\bar{j})){s_{A,n}}}{s_{A,n}}^{\bar{i}+\alpha(\bar{k}-1-\bar{j})}A^{\bar{k}-1+2\theta(\bar{k}-1-\bar{j})}$
$\displaystyle\times\left(\sum_{(k,j,i)\in\bar{E}}\left(1+\frac{\tau}{{s_{A,n}}}\right)^{i}(k-j)\tilde{C}_{k,j,i}L^{k-1-j}+o(1)\right).$

Now, for any $\delta_{0}>0$, we consider $\tau_{n}(\delta_{0})>0$ such that

$e^{\frac{\tau}{2}}e^{(1-\frac{\bar{k}}{2}-\theta(\bar{k}-1-\bar{j})){s_{A,n}}}{s_{A,n}}^{\bar{i}+\alpha(\bar{k}-1-\bar{j})}A^{\bar{k}-1+2\theta(\bar{k}-1-\bar{j})}=\delta_{0}.$

Clearly, it holds that

$\tau_{n}(\delta_{0})\sim[\bar{k}-2+\theta(\bar{k}-1-\bar{j})]{s_{A,n}}\mbox{
as }n\to\infty.$

Therefore, we see that

$v_{b,1,0}({s_{A,n}}+\tau_{n}(\delta_{0}))=\delta_{0}\left(\sum_{(k,j,i)\in\bar{E}}\left([\bar{k}-1+\theta(\bar{k}-1-\bar{j})\right)^{i}(k-j)\tilde{C}_{k,j,i}L^{k-1-j}+o(1)\right).$

This forces the coefficient of $\delta_{0}$ to be zero:

$\sum_{(k,j,i)\in\bar{E}}\left([\bar{k}-1+\theta(\bar{k}-1-\bar{j})\right)^{i}(k-j)\tilde{C}_{k,j,i}L^{k-1-j}=0,$

otherwise, $v_{b,1,0}({s_{A,n}}+\tau_{n}(\delta_{0}))$ will be large as in
(4.21), which leads to a contradiction.
Since $L\neq 0$ and $\tilde{C}_{k,j,i}\neq 0$, for all $(k,j,i)\in\bar{E}$
(remember that $\bar{E}\subset H_{2}$ defined in (5.17)), this identity
contains at least two terms, and this is precisely the desired polynomial
relation.
This concludes the proof of Proposition 5.4 and finishes Part 1 dedicated to
the super-quadratic case. ∎

Part 2: The quadratic regime

In this part, we consider the quadratic regime, where

${a_{n,1}}\sim L{a_{n,2}}^{2}\mbox{ as }n\to\infty,$ (5.35)

for some $L>0$, and show that $L$ satisfies a polynomial equation whose
coefficients are given by the Taylor expansion of the solution.

Note first from the notation (4.11) and the definition (4.9) of $B_{n}$ that
we have

$B_{n}\sim LAe^{-\frac{{s_{A,n}}}{2}}\mbox{ as }n\to\infty.$

The proof is in fact a simple adaptation of our argument in the super-
quadratic regimes, given in Part 1. Let us then follow that part step by step,
and see what changes.

In Part 1, the outcome of Step 1 is the following:

$C_{m+1,m}+LC_{m,m-2}=0.$ (5.36)

If $C_{m,m-2}\neq 0$, then, we have our polynomial and we are done (in fact,
we have more, in the sense that $L$ enjoys only one value :
$-C_{m+1,m}/C_{m,m-2}$).
If $C_{m,m-2}=0$, then we also have $C_{m,m-3}=0$, because the multilinear
form in (1.10) is nonpositive. In other words, we have exactly the same
conclusions as in Step 1 of Part 1 in the super-quadratic case (in fact, we
have more, since $C_{m+1,m}=0$ from (5.36)). For short, we can carry on all
the next steps of the super-quadratic case up to to the end of Part 1, and see
that $L$ satisfies a polynomial equations. For the reader’s convenience, we
would like to mention that hypothesis (5.35) makes many steps either non
relevant or trivial:
\- In Step 2 of Part 1, $\theta=\frac{1}{2}$ which makes Proposition 5.1 non
relevant. Moreover, (5.35) is stronger than (5.5).
\- In Step 3, estimate (5.35) is stronger than (5.20), and $\alpha$ defined in
(5.21) is zero. Accordingly, Proposition 5.3 is non relevant.
\- As for Step 4, again, estimate (5.35) is stronger than (5.27), and the
first assertion of Proposition 5.4 is clear ($\psi_{n}\to L$ as $\to\infty$).
\- In conclusion, only the second assertion of Proposition 5.4 remains
relevant, and provides us with the polynomial relation for $L$.

Part 3: Conclusion of the proof of Theorem 3

From the study of the super-quadratic case in Part 1, and also the quadratic
case in Part 2, we see that either we are in the subquadratic case
${a_{n,1}}=o({a_{n,2}})^{2}$, or

$B_{n}\sim LA^{2\theta}e^{-\theta{s_{A,n}}}{s_{A,n}}^{\alpha}\mbox{ as
}n\to\infty,$

where $\theta$ and $\alpha$ enjoy only a finite set of rational values, and
$L>0$ is a solution of a polynomial equations whose coefficients depend on the
Taylor expansion of the solution. By definitions (4.11) and (4.9) of
${s_{A,n}}$ and $B_{n}$, we see that

${a_{n,1}}\sim L{a_{n,2}}^{2\theta+1}|2\log{a_{n,2}}|^{\alpha},$

with $2\theta+1\in(0,1]$, which is the desired estimate in Theorem 3 The set
where $2\theta+1$ lives directly follows from item (iii) in Proposition 5.1.

This concludes the proof of Theorem 3. ∎

## References

  * [1] X. Y. Chen and H. Matano.  Convergence, asymptotic periodicity, and finite-point blow-up in one-dimensional semilinear heat equations.  J. Differential Equations, 78(1):160–190, 1989.
  * [2] S. Filippas and R. V. Kohn.  Refined asymptotics for the blowup of $u_{t}-{\Delta}u=u^{p}$.  Comm. Pure Appl. Math., 45(7):821–869, 1992.
  * [3] S. Filippas and W. X. Liu.  On the blowup of multidimensional semilinear heat equations.  Ann. Inst. H. Poincaré Anal. Non Linéaire, 10(3):313–344, 1993\.
  * [4] T. Ghoul, V.T. Nguyen, and H. Zaag.  Refined regularity of the blow-up set linked to refined asymptotic behavior for the semilinear heat equation.  Adv. Nonlinear Stud., 17(1):31–54, 2017.
  * [5] Y. Giga and R. V. Kohn.  Asymptotically self-similar blow-up of semilinear heat equations.  Comm. Pure Appl. Math., 38(3):297–319, 1985.
  * [6] Y. Giga and R. V. Kohn.  Nondegeneracy of blowup for semilinear heat equations.  Comm. Pure Appl. Math., 42(6):845–884, 1989.
  * [7] Y. Giga, S. Matsui, and S. Sasayama.  Blow up rate for semilinear heat equations with subcritical nonlinearity.  Indiana Univ. Math. J., 53(2):483–514, 2004.
  * [8] M. A. Herrero and J. J. L. Velázquez.  Blow-up profiles in one-dimensional, semilinear parabolic problems.  Comm. Partial Differential Equations, 17(1-2):205–219, 1992.
  * [9] M. A. Herrero and J. J. L. Velázquez.  Blow-up behaviour of one-dimensional semilinear parabolic equations.  Ann. Inst. H. Poincaré Anal. Non Linéaire, 10(2):131–189, 1993\.
  * [10] S. Khenissy, Y. Rébaï, and H. Zaag.  Continuity of the blow-up profile with respect to initial data and to the blow-up point for a semilinear heat equation.  Ann. Inst. H. Poincaré Anal. Non Linéaire, 28(1):1–26, 2011\.
  * [11] F. Merle, P. Raphaël, and J. Szeftel.  On strongly anisotropic type I blowup.  Int. Math. Res. Not., 2020(2):541–606, 2020.
  * [12] F. Merle and H. Zaag.  Optimal estimates for blowup rate and behavior for nonlinear heat equations.  Comm. Pure Appl. Math., 51(2):139–196, 1998.
  * [13] F. Merle and H. Zaag.  A Liouville theorem for vector-valued nonlinear heat equations and applications.  Math. Annalen, 316(1):103–137, 2000.
  * [14] R. Schweyer.  Type II blow-up for the four dimensional energy critical semi linear heat equation.  J. Funct. Anal., 263(12):3922–3983, 2012.
  * [15] J. J. L. Velázquez.  Higher-dimensional blow up for semilinear parabolic equations.  Comm. Partial Differential Equations, 17(9-10):1567–1596, 1992\.
  * [16] J. J. L. Velázquez.  Classification of singularities for blowing up solutions in higher dimensions.  Trans. Amer. Math. Soc., 338(1):441–464, 1993.
  * [17] H. Zaag.  On the regularity of the blow-up set for semilinear heat equations.  Ann. Inst. H. Poincaré Anal. Non Linéaire, 19(5):505–542, 2002\.
  * [18] H. Zaag.  One dimensional behavior of singular ${N}$ dimensional solutions of semilinear heat equations.  Comm. Math. Phys., 225(3):523–549, 2002.
  * [19] H. Zaag.  Regularity of the blow-up set and singular behavior for semilinear heat equations.  In Mathematics & mathematics education (Bethlehem, 2000), pages 337–347. World Sci. Publishing, River Edge, NJ, 2002.
  * [20] H. Zaag.  Determination of the curvature of the blow-up set and refined singular behavior for a semilinear heat equation.  Duke Math. J., 133(3):499–525, 2006.

Address:
CY Cergy Paris Université, Département de mathématiques, 2 avenue Adolphe
Chauvin, BP 222, 95302 Cergy Pontoise cedex, France.

    e-mail<EMAIL_ADDRESS>

Université Sorbonne Paris Nord, Institut Galilée, Laboratoire Analyse,
Géométrie et Applications, CNRS UMR 7539, 99 avenue J.B. Clément, 93430
Villetaneuse, France.

    e-mail<EMAIL_ADDRESS>


# Ensemble-based gradient inference for particle methods in optimization and
sampling

Claudia Schillings, Claudia Totzeck and Philipp Wacker

###### Abstract

We propose an approach based on function evaluations and Bayesian inference to
extract higher-order differential information of objective functions from a
given ensemble of particles. Pointwise evaluation $\\{V(x^{i})\\}_{i}$ of some
potential $V$ in an ensemble $\\{x^{i}\\}_{i}$ contains implicit information
about first or higher order derivatives, which can be made explicit with
little computational effort (ensemble-based gradient inference – EGI). We
suggest to use this information for the improvement of established ensemble-
based numerical methods for optimization and sampling such as Consensus-based
optimization and Langevin-based samplers. Numerical studies indicate that the
augmented algorithms are often superior to their gradient-free variants, in
particular the augmented methods help the ensembles to escape their initial
domain, to explore multimodal, non-Gaussian settings and to speed up the
collapse at the end of optimization dynamics. The code for the numerical
examples in this manuscript can be found in the paper’s Github repository111
https://github.com/MercuryBench/ensemble-based-gradient.git

##  1 Introduction

Global optimization of nonconvex objective functions and sampling from
nonstandard distributions are widespread applications in industry, finance and
other disciplines. Although the problems read very simple, and there is a
variety of algorithms available, it is still very challenging to design well-
performing algorithms and to prove their convergence. We motivate both the
optimzation and sampling problem via the inverse setting, where the goal is to
identify unknown parameters $x\in\mathbb{R}^{d}$ from noisy observations

$y=G(x)+\eta\,,$ (1.1)

with $G:\mathbb{R}^{d}\to\mathbb{R}^{k}$ denoting the parameter-to-observation
map and $\eta$ being the noise. Under suitable assumptions on $G$ and $\eta$,
application of Bayes’ theorem gives as solution of the (Bayesian) inverse
problem the following characterization of the posterior distribution

$\frac{\mathrm{d}\mu(x)}{\mathrm{d}\lambda(x)}\propto\exp(-V(x))$ (1.2)

w.r. to the Lebesgue measure $\lambda$, where $V$ denotes a regularized
potential. The aim of sampling methods is then to generate samples
(approximately) distributed according to the posterior distribution. The
connection to the optimization setting is via the maximum-a-posteriori (MAP)
estimate (assuming the unique existence) given by

$x^{\star}=\operatorname{arg\,max}\limits_{x}\exp(-V(x))=\operatorname{arg\,min}\limits_{x}V(x)\,.$
(1.3)

The Bayesian approach to inverse problems has become very popular over the
last decades and there has been a lot of research effort towards efficient
methods, in particular gradient-free methods in order to allow the use of
black-box solvers for the underlying forward problem. Kalman-Wasserstein flows
have been the basis to design various efficient, gradient-free samplers, e.g.
Ensemble Kalman Sampler Garbuno-Inigo et al. (2020a); Nüsken and Reich (2019),
Affine Invariant Interacting Langevin Dynamics Garbuno-Inigo et al. (2020b).
Also in the optimization setting, gradient free methods such as Consensus-
based optimization method (CBO) and Ensemble Kalman Inversion (EKI) have
become very popular, see e.g. Carrillo et al. (2018, 2021); Totzeck (2022);
Schillings and Stuart (2017); Blömker et al. (2019) and the references
therein. These methods have in common that they rely on an ensemble of
particles, which is transformed into posterior samples or concentrates around
the minimizer in the limit ’time to infinity’. We will show in the following
that the ensemble itself can be used to efficiently approximate derivative
information (without any further evaluations of the forward problem). In order
to do so, we will formulate the approximation of derivative information as an
inverse problem and extract the information from given function evaluations.
To demonstrate the potential of the simultaneous estimation of the derivative
information and parameter estimation, we consider the ensemble based gradient
augmentation of the following state-of-the art optimization and sampling
methods:

  * •

Consensus-based optimization (CBO)

  * •

Ensemble Langevin sampler (LS)

  * •

Metropolis-adjusted Langevin algorithm (MALA)

  * •

Affine invariant interacting Langevin dynamics for Bayesian inference (ALDI)
and Ensemble Kalman sampler (EKS).

In order to keep the presentation self-contained and for later reference we
review the methods in the following.

###  1.1 Review of existing computational methods

We consider a potential $V$ as described above which may be the negative
logdensity of some Bayesian posterior measure, or just a function we want to
minimize. We will call the task of finding the minimal value of $V$ (if it
exists) the optimization task, and the task of constructing samples from the
measure with unnormalized Lebesgue-density $\exp(-V)$ will be called the
sampling task.

#### CBO (for optimization)

The consensus-based optimization (CBO) method Carrillo et al. (2018); Pinnau
et al. (2017) describes the collective dynamical behavior of an ensemble of
particles exploring the state space experiencing some diffusive behaviour and
finally collapsing onto a joint weighted ensemble mean, which is by
construction an approximation of the minimizer of $V.$ Given a measure
$\rho\in\mathcal{M}(\mathbb{R}^{d})$ we define the weighted mean of $\rho$ as

$\displaystyle\mathsf{m}_{\beta}(\rho):=\frac{\int x\exp(-\beta
V(x))\,\mathrm{d}\rho(x)}{\int\exp(-\beta V(x))\,\mathrm{d}\rho(x)}.$ (1.4)

This definition of weighted mean is motivated by the Laplace principle Pinnau
et al. (2017). Indeed, assuming that $V$ admits a unique global minimizer, it
can be shown that for $\beta\rightarrow\infty,$ $\mathsf{m}_{\beta}(\rho)$
tends to the
$\operatorname{arg\,min}\limits_{x\in\operatorname{supp}{\rho}}V(x)$.

The CBO dynamics combines exploration of the landscape given by $V$ with
aggregation of the ensemble around $\mathsf{m}_{\beta}(\rho).$ More
succinctly, it follows the following system of stochastic differential
equations (SDEs), where $\rho:=\sum_{i=1}^{J}\delta_{x^{i}}.$

$\displaystyle\,\mathrm{d}x^{i}$
$\displaystyle=-\lambda(x^{i}-\mathsf{m}_{\beta}(\rho))\,\mathrm{d}t+\sigma\left|x^{i}-\mathsf{m}_{\beta}(\rho)\right|\,\mathrm{d}W_{t}^{i},\qquad
i=1\ldots,J$ (1.5)

with initial conditions $x_{0}^{i}\sim\mathcal{P}_{2}(\mathbb{R}^{d})$ drawn
independently. It is important to note that both the drift and the diffusion
part of the dynamics scale with the distance
$|x^{i}-\mathsf{m}_{\beta}(\rho)|.$ On the one hand this allows for collapse
at $\mathsf{m}_{\beta}(\rho)$ which is important for optimization tasks. On
the other hand, the dynamics slows down as the ensemble variance becomes
small.

We will call the algorithm corresponding to (1.5) “vanilla CBO” and in case
the componentwise noise modification of Carrillo et al. (2018) is used, we
will specify this as “component-wise noise vanilla CBO”. In Pinnau et al.
(2017) the Consensus-based optimization method (CBO) was proposed as
alternative to heuristic gradient-free particle optimization methods such as
evolutionary or genetic algorithms and simulated annealing (Simon, 2013; Back,
1996; van Laarhoven and Aarts, 1987) with the intention to prove convergence
for the corresponding mean-field dynamics (Carrillo et al., 2021). Based on
these publications a variety of researchers from stochastic analysis (Huang
and Qiu, 2022; Kalise et al., 2022), PDE analysis (Fornasier et al., 2021b,
2022) found ways to relax assumptions, give further insight in the internal
process of the method (Fornasier et al., 2021b), open the class of
applications like constraint problems (Fornasier et al., 2021a, 2020; Ha et
al., 2022), machine learning (Carrillo et al., 2021), multiobjective problems
(Borghi et al., 2022; Klamroth et al., 2022) and sampling (Carrillo et al.,
2022b). Others aim to make the algorithm more reliable (Totzeck and Wolfram,
2020). As the original method with minimal dynamics is fairly understood,
there are mainly two paths to go, either improve the method itself or use the
established ideas to prove convergence results for other well-known
optimization dynamics, see for example (Grassi and Pareschi, 2021; Huang,
2021) for Particle Swarm Optimization.

#### Langevin dynamics and MALA (for sampling)

A well-known sampling dynamics is given by the (overdamped) Langevin equation,
cp. Yang et al. (2019) for a recent overview for the connection of Langevin
dynamics and Markov chain Monte Carlo (MCMC) methods. Given a symmetrical and
positive definite matrix $M$, the (overdamped) Langevin equation
preconditioned by $M$ models a particle subject to the influence of the
potential $V$:

$\,\mathrm{d}x_{t}=-M\nabla
V(x_{t})\,\mathrm{d}t+\sqrt{2}M^{1/2}\,\mathrm{d}W_{t}.$ (1.6)

It can be shown that the invariant measure of this stochastic differential
equation is equal to $\mu$, i.e. the target measure $\mu$ with unnormalized
Lebesgue-density $\exp(-V)$ Garbuno-Inigo et al. (2020b); Jordan et al.
(1998). In applications the gradient of $V$ may be expensive to compute or
simply not available. It is therefore of interest to construct a surrogate
model with a gradient approximation.

#### ALDI and EKS

Another recent ensemble-based method with pointwise function evaluations is
the gradient-free ALDI sampler Garbuno-Inigo et al. (2020b). It performs an
implicit approximation to the gradient of the loglikelihood. There is a close
relationship between ALDI and the Ensemble Kalman sampler (EKS) Garbuno-Inigo
et al. (2020a); Nüsken and Reich (2019) which is an efficient gradient-free
interacting particle sampler originating from the Ensemble Kalman filter. EKS
can be applied if the negative logdensity of $\rho$ can be written in the form
$V(x)=\frac{1}{2}\|y-A(x)\|_{\Gamma}^{2}+\frac{1}{2}\|x-\mu_{0}\|_{\Sigma_{0}}^{2}$,
which is typical for a Bayesian statistical problem. In this special case, the
gradient-free ALDI is the coupled system of SDEs given by

$\begin{split}\,\mathrm{d}x^{i}_{t}&=-\\{D(X_{t})\Gamma^{-1}(A(x^{i}_{t})-y)+C(X_{t})\Sigma_{0}^{-1}(x_{t}^{i}-\mu_{0})\\}\,\mathrm{d}t\\\
&\qquad+\frac{d+1}{J}(x_{t}^{i}-\bar{x}_{t})\,\mathrm{d}t+\sqrt{2}C(X_{t})^{\frac{1}{2}}\,\mathrm{d}W_{t}^{i}\end{split}$
(1.7)

where $d$ is the dimension of the space and

$\begin{split}C(X)&=\frac{1}{N}\sum_{j=1}^{N}(x^{j}-\bar{x})\otimes(x^{j}-\bar{x})\\\
D(X)&=\frac{1}{N}\sum_{j=1}^{N}(x^{j}-\bar{x})\otimes(A(x^{j})-\bar{A}),\quad\bar{A}=\frac{1}{N}\sum_{j=1}^{N}A(x^{j})\end{split}$
(1.8)

with initial conditions analogous to the methods above. If $A$ is a linear
map, samples from $\rho$ are provable accurately generated. One of its main
strengths stems from the fact that the Ensemble Kalman method constructs a
covariance-preconditioned approximation to the gradient: Note that if
$V(x)=\frac{1}{2}\|A(x)-y\|_{\Gamma}^{2}$ and $A$ linear, then

$\displaystyle C(X)\cdot\nabla V(x^{i})$
$\displaystyle=\frac{1}{N}\sum_{j=1}^{N}(x^{j}-\bar{x})\langle
x^{j}-\bar{x},\nabla V(x^{i})\rangle$
$\displaystyle=\frac{1}{N}\sum_{j=1}^{N}(x^{j}-\bar{x})\langle
x^{j}-\bar{x},A^{T}\Gamma^{-1}(Ax^{i}-y)\rangle$
$\displaystyle=\frac{1}{N}\sum_{j=1}^{N}(x^{j}-\bar{x})\langle
A(x^{j})-A(\bar{x}),\Gamma^{-1}(Ax^{i}-y)\rangle$
$\displaystyle=D(X)\cdot\Gamma^{-1}(A(x^{i})-y).$

For $A$ linear, the approximation is exact. In fact, in the linear regime
(1.7) corresponds to

$\,\mathrm{d}x_{t}^{i}=-C(X)\cdot\nabla
V(x_{t}^{i})\,\mathrm{d}t+\frac{d+1}{J}(x_{t}^{i}-\bar{x}_{t})\,\mathrm{d}t+\sqrt{2}C(X_{t})^{\frac{1}{2}}\,\mathrm{d}W_{t}^{i},$

which is an overdamped Langevin equation with a correction term ensuring the
affine invariance of the resulting scheme Garbuno-Inigo et al. (2020b).

###  1.2 Main idea and contributions of the paper

We argue that tracking an ensemble of particles $\\{x^{i}\\}_{i=1}^{J}$ with
pointwise evaluations $\\{V(x^{i})\\}_{i=1}^{J}$ of some potential function
$V$ carries implicit gradient (and higher-order differential) information
which can be converted into explicit information basically “for free”, via an
often negligible cost of solving one linear equation system. This can and
should be leveraged in various optimization and sampling tasks:

The CBO algorithm for optimization is completely gradient-agnostic and can
benefit from incorporation of gradient information via EGI in order to speed
up local convergence (“better convergence to minima”) and improve overall
performance (“convergence to better minima”).

Langevin dynamics and MALA in its original form need explicit gradient
information. By managing a full ensemble instead of just one particle and
thereby providing approximated gradient information for each member of the
ensemble via EGI, we can run approximated Langevin dynamics and MALA without
the need for an explicit gradient functional.

EKS and ALDI utilize a powerful property of the empirical covariance of an
ensemble in order to compute a preconditioned gradient without the need for
any additional computation. This allows us to sample very efficiently from a
Bayesian posterior. Unfortunately, this is an approximation that holds only in
the linear, Gaussian setting. With EGI, we can obtain a better approximation
of the gradient, allowing for construction of non-Gaussian posterior samples.

Summarized, we present, motivate, and test the following novel algorithms.

  * •

EGI: Gradient approximation from ensemble point evaluation (as a general
method)

  * •

EGI-CBO: CBO augmented by gradient information

  * •

EGI-LS: Coupled Langevin dynamics with approximated gradient information

  * •

EGI-MALA: Coupled MALA dynamics with approximated gradient information

  * •

EGI-ALDI: Gradient-free ALDI with approximated gradients, with some variants

The rest of the article is structured as follows: Section 2 describes how
pointwise evaluation of some function $V$ evaluated in an ensemble
$\\{x^{i}\\}_{i=1}^{J}$ can be converted into explicit inexact gradient
information via solving a linear inverse problem from straightforward Taylor
approximation. We call this method EGI – Ensemble-based gradient inference.
Section 3 develops EGI-CBO, a consensus-based optimization method augmented by
inexact gradient information. Finally, Section 4 generalizes this idea to
ensemble-based sampling methods which normally use exact gradient evaluation,
but can be shown to work with inexact gradients supplied by EGI as well.

##  2 Ensemble-based gradient inference (EGI) from pointwise ensemble
evaluation

Consider an unknown function $V:\mathbb{R}^{d}\to\mathbb{R}$ evaluated
pointwise in an ensemble $\\{x^{k}\\}_{k=1}^{J}\subset\mathbb{R}^{d}$. Our
goal is to infer gradient information (as well as possibly higher-order
derivatives) at one or each of the ensemble members’ positions $x^{j}$, i.e.
$\nabla V(x^{j})$ (and possibly higher-order derivatives). In short, we want
to solve the ill-posed problem

$\\{V(x^{k})\\}_{k=1}^{J}\mapsto\nabla V(x^{\star}),$

where $x^{\star}$ is an element of the ensemble $\\{x^{k}\\}_{k=1}^{J}$222It
is conceivable to allow $x^{\star}$ to be an arbitrary point “near” the
ensemble by appropriately modifying the ansatz below.. We now describe our
approach where we set $x^{\star}=x^{j}$ for some fixed $j\in\\{1,\ldots,J\\}$.
By Taylor’s formula, for any $i\in\\{1,\ldots,J\\}$,

$\displaystyle V(x^{i})-V(x^{j})$ $\displaystyle=\nabla
V(x^{j})^{T}(x^{i}-x^{j})$
$\displaystyle+\frac{1}{2}(x^{i}-x^{j})^{T}HV(x^{j})(x^{i}-x^{j})+\varepsilon_{i}^{j}\cdot\frac{\|x^{i}-x^{j}\|^{3}}{6},$

where we model the action of the unknown tensors $D^{3}V(x^{j})$ on
$[x^{i}-x^{j},x^{i}-x^{j},x^{i}-x^{j}]$ by independent random variables
$\\{\varepsilon_{i}^{j}\\}_{i=1}^{J}\sim N(0,\gamma^{2})$.We make the ansatz
$\nabla
V(x^{j})\approx\sum_{k=1}^{J}u^{j,1}_{k}\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}$
and
$HV(x^{j})\approx\sum_{k=1}^{J}u_{k}^{j,2}\frac{(x^{k}-x^{j})}{\|x^{k}-x^{j}\|}\otimes\frac{(x^{k}-x^{j})}{\|x^{k}-x^{j}\|}$
for $u^{j,1}_{k},u^{j,2}_{k}\in\mathbb{R}$ yet to be determined. Clearly, this
means that our approximation for $\nabla$ will be an element of the affine
space spanned by the ensemble, also $HV$ is now approximated by a rank-$J$
matrix. We introduce an additional error slack term $\xi>0$, which will be
used for localization purposes, see Remark 2.3. Then

$\displaystyle V(x^{i})-V(x^{j})$
$\displaystyle=\sum_{k=1}^{J}u_{k}^{j,1}\cdot\frac{(x^{k}-x^{j})^{T}(x^{i}-x^{j})}{\|x^{k}-x^{j}\|}$
$\displaystyle+\frac{1}{2}\sum_{k=1}^{J}u_{k}^{j,2}\cdot\left[\frac{(x^{k}-x^{j})^{T}(x^{i}-x^{j})}{\|x^{k}-x^{j}\|}\right]^{2}+\varepsilon_{i}^{j}\cdot\left(\frac{\|x^{i}-x^{j}\|^{3}}{6}+\xi\right).$

######  Remark 2.1.

For consistency, we set $0/0=0$ here.

######  Remark 2.2.

Alternatively we could include mixed terms in the formulation of $HV(x^{j})$,
i.e.
$HV(x^{i})\approx\sum_{k,l}u_{k,l}^{j,2}\frac{(x^{k}-x^{j})}{\|x^{k}-x^{j}\|}\otimes\frac{(x^{l}-x^{j})}{\|x^{l}-x^{j}\|}$
for $u_{k,l}^{j,2}\in\mathbb{R}.$ Note that this comes at the cost that the
number of coefficients scales with $J^{2}$ instead of $J$. For simplicity of
presentation, we stick with the rank-$J$ ansatz for $HV(x^{j})$.

######  Remark 2.3.

The motivation for the introduction of the error slack term $\xi$ is two-fold:
First, if there are two ensemble members $x^{i},x^{j}$ very close to each
other, then this error term acts as a safeguard against unwanted overfitting.
Second, if all ensemble members are very close to each other (in the magnitude
of machine precision), then it buffers unwanted numerical instability. This
error term can also be removed from the model by setting $\xi=0$. On the other
hand, letting $\xi\to\infty$ recovers the least-squares quadratic regression
function through the data points, under the additional condition that the
regression function passes through $(x^{j},V(x^{j}))$. In a sense, $\xi=0$
corresponds to “local” gradient information, because it takes into account
mostly nearby points, and $\xi\to\infty$ yields “global” gradient information.
The effect of $\xi$ is illustrated in Figure 1.

We now describe how to set up gradient inference as a linear inverse problem.

###  2.1 Ensemble-based gradient inference as an inverse problem

Let $\\{x^{i}\\}_{i=1}^{J}$ be an ensemble of points with pointwise
evaluations $\\{V(x^{i})\\}_{i=1}^{J}$ of an unknown function $V$. We fix
$j\in\\{1,\ldots,J\\}$, and our goal is to infer $\nabla V(x^{j})$ (and
$HV(x^{j})$) from the data. We set the following notational shorthands:

$\displaystyle X^{j}$
$\displaystyle=(x^{1}-x^{j},\ldots,x^{i}-x^{j},\ldots,x^{J}-x^{j})\in\mathbb{R}^{d\times
J},\qquad$ $\displaystyle X^{j}_{:,i}=x^{i}-x^{j},$ $\displaystyle Z^{j}$
$\displaystyle=\left(\frac{x^{1}-x^{j}}{\|x^{1}-x^{j}\|},\ldots,\frac{x^{J}-x^{j}}{\|x^{J}-x^{j}\|}\right)\in\mathbb{R}^{d\times
J},\qquad$ $\displaystyle Z_{:,i}^{j}=\frac{x^{i}-x^{j}}{\|x^{i}-x^{j}\|},$
$\displaystyle y^{j}$
$\displaystyle=(V(x^{1})-V(x^{j}),\ldots,V(x^{J})-V(x^{j}))^{T}\in\mathbb{R}^{J},\qquad$
$\displaystyle y_{i}^{j}=V(x^{i})-V(x^{j}).$

Then with $u^{j,i}=(u_{k}^{j,i})_{k=1}^{J}\in\mathbb{R}^{J}$ and
$u^{j}=(u^{j,1},u^{j,2})\in\mathbb{R}^{2J}$,

$\displaystyle y^{j}$
$\displaystyle=\begin{bmatrix}{X^{j}}^{T}Z^{j}&\frac{({X^{j}}^{T}Z^{j})^{\odot
2}}{2}\end{bmatrix}\cdot\begin{pmatrix}u^{j,1}\\\
u^{j,2}\end{pmatrix}+\left(\frac{1}{3!}\mathrm{diag}(\\{\|x^{i}-x^{j}\|^{3}\\}_{i=1}^{J})+\xi\right)\cdot\varepsilon^{j},$
(2.1)

where the Hadamard square $({X^{j}}^{T}Z^{j})^{\odot 2}$ is taken component-
wise, i.e. $[({X^{j}}^{T}Z^{j})^{\odot
2}]_{i,j}=[({X^{j}}^{T}Z^{j})_{i,j}]^{2}$. We can further compact this form by
setting $A^{j}=({X^{j}}^{T}Z^{j},\frac{({X^{j}}^{T}Z^{j})^{\odot
2}}{2})\in\mathbb{R}^{d\times(2J)}$ and
$\Gamma_{\xi,\gamma}^{j}=\left(\frac{\gamma^{2}}{3!}\mathrm{diag}(\\{\|x^{i}-x^{j}\|^{3}\\}_{i=1}^{J})+\xi\gamma^{2}\right)$.
Then

$y^{j}=A^{j}u^{j}+\varepsilon^{j}$ (2.2)

with $\varepsilon^{j}\sim N(0,\Gamma_{\xi,\gamma}^{j})$ is a linear inverse
problem for the coefficient vector $u^{j}$. If ${\bar{u}}^{j}$ is the least
squares solution of this inverse problem, then we can recover gradient
information via

$\displaystyle G^{j}_{\xi}(\\{x^{i}\\}_{i=1}^{J})$
$\displaystyle:=Z^{j}\cdot{{\bar{u}}^{j,1}}=\sum_{k=1}^{J}{\bar{u}}_{k}^{j,1}\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}\approx\nabla
V(x^{j})$ (2.3) $\displaystyle H^{j}_{\xi}(\\{x^{i}\\}_{i=1}^{J})$
$\displaystyle:=\sum_{k=1}^{J}{\bar{u}}_{k}^{j,2}\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}\otimes\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}\approx
HV(x^{j})$ (2.4)

This notation hides the implicit dependence on $V$, and the specific way
$u^{j}$ is obtained (although we assume this to be the least square solution
to (2.2) unless stated otherwise). For brevity we write
$G^{j}(\\{x^{i}\\}_{i=1}^{J}):=G^{j}_{0}(\\{x^{i}\\}_{i=1}^{J})$ in the case
$\xi=0$ and similar for the Hessian $H^{j}$. Details of the inference are
given in Algorithm 1. In the following we will refer to the case $\xi=0$ as
local approximation and $\xi>0$ global gradient approximation. This
characterization is motivated by the observation that $\xi$ enters the error
term in (2.1). Therefore, for $\xi=0$, the Taylor approximation in the
location of a particle $x^{i}$ close to $x^{j}$ needs to be much better than
in any particle far away from $x^{j}$. In the other extreme case
$\xi\rightarrow\infty$ all particles enter with equal weight in the gradient
approximation, hence the second-order polynomial approximation will be a
“global” approximation instead of a local Taylor approximation. See Figure 1
for a visualization of the effect of $\xi$.

######  Remark 2.4.

The Hessian approximation can be written compactly in Python Einstein sum
notation via $H^{j}(\\{x^{i}\\}_{i=1}^{J})=$ np.einsum(’k,lk,mk->lm’,
${\bar{u}}^{j,2}$, $Z^{j}$, $Z^{j}$).

Data: $\\{x^{i},V^{i}\\}_{i=1}^{J},\gamma>0,\xi\geq 0$, $j\in\mathbb{N}$

Result:
$G=G^{j}_{\xi}(\\{x^{i}\\}_{i=1}^{J}),H=H^{j}_{\xi}(\\{x^{i}\\}_{i=1}^{J})$
(approx. to $\nabla V(x^{j})$ and $HV(x^{j})$)

1 $X\leftarrow(x^{1}-x^{j},\ldots,x^{J}-x^{j})\in\mathbb{R}^{d\times J}$

2
$Z\leftarrow\begin{pmatrix}\frac{x^{1}-x^{j}}{\|x^{1}-x^{j}\|},&\cdots&,\frac{x^{J}-x^{j}}{\|x^{J}-x^{j}\|}\end{pmatrix}\in\mathbb{R}^{d\times
J}$

3 $y\leftarrow(V^{1}-V^{j},\ldots,V^{J}-V^{j})\in\mathbb{R}^{J}$

4 $A\leftarrow(X^{T}Z,\tfrac{1}{2}(X^{T}Z)^{\odot
2})\in\mathbb{R}^{d\times(2J)}$

5
$\Gamma\leftarrow\gamma^{2}\cdot\left[\tfrac{1}{3!}\operatorname{diag}(\|x^{i}-x^{j}\|^{3})_{i=1}^{J}+\xi\right]\in\mathbb{R}^{d\times
d}$

$(u^{1},u^{2})^{T}\leftarrow\mathtt{lsq}(\Gamma^{-1}A,\Gamma^{-1}y)$// Solve
$\Gamma^{-1}A(u^{1},u^{2})^{T}=\Gamma^{-1}y$ via Least-Squares

6 $G\leftarrow\sum_{k=1}^{J}{u}_{k}^{1}\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}$

7
$H\leftarrow\sum_{k=1}^{J}{u}_{k}^{2}\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}\otimes\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}$

Algorithm 1 Ensemble-based Gradient Inference (EGI)

######  Remark 2.5.

The least-squares approach can without major changes be turned into a Bayesian
sampling algorithm for gradient and Hessian information by solving the linear
inverse problem (2.2) via the Gaussian update formula, taking samples
$\\{u_{n}^{j}\\}_{n}$, and propagating the samples $u^{j}_{n}$ via a mapping
similar to Equations 2.3 and 2.4. This is described in Algorithm 2, but we
forgo fixing any notation for this in order to keep the presentation clearer.
We just remark that any subsequent algorithm can be modified to use Bayesian
samples of the approximated gradient (Algorithm 2) instead of the
deterministic least squares approximation Algorithm 1, as will be done in
Algorithm 8.

Data: $\\{x^{i},V^{i}\\}_{i=1}^{J},\gamma>0,\xi\geq 0$, index
$j\in\mathbb{N}$, $\Sigma>0$, number of samples $N\in\mathbb{N}$

Result: samples $\\{G^{n,j}\\}_{n=1}^{N}$, $\\{H^{n,j}\\}_{n=1}^{N}$
(approximations to $\nabla V(x^{j})$ and $HV(x^{j})$) and Bayesian maximum-a-
posteriori estimators $G_{\text{MAP}},H_{\text{MAP}}$.

1 $X\leftarrow(x^{1}-x^{j},\ldots,x^{J}-x^{j})\in\mathbb{R}^{d\times J}$

2
$Z\leftarrow\begin{pmatrix}\frac{x^{1}-x^{j}}{\|x^{1}-x^{j}\|},&\cdots&,\frac{x^{J}-x^{j}}{\|x^{J}-x^{j}\|}\end{pmatrix}\in\mathbb{R}^{d\times
J}$

3 $y\leftarrow(V^{1}-V^{j},\ldots,V^{J}-V^{j})\in\mathbb{R}^{J}$

4 $A\leftarrow(X^{T}Z,\tfrac{1}{2}(X^{T}Z)^{\odot
2})\in\mathbb{R}^{d\times(2J)}$

5
$\Gamma\leftarrow\gamma^{2}\cdot\left[\tfrac{1}{3!}\operatorname{diag}(\|x^{i}-x^{j}\|^{3})_{i=1}^{J}+\xi\right]\in\mathbb{R}^{d\times
d}$

6 $K\leftarrow\Sigma A^{T}(\Gamma+A\Sigma A^{T})^{-1}$

7 $\hat{\mu}\leftarrow(\hat{\mu}^{1},\hat{\mu}^{2})\leftarrow Ky$

8 $\hat{\Sigma}\leftarrow\Sigma-KA\Sigma$

9
$G_{\text{MAP}}\leftarrow\sum_{k=1}^{J}{\hat{\mu}}^{1}_{k}\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}$

10
$H_{\text{MAP}}\leftarrow\sum_{k=1}^{J}{\hat{\mu}}^{2}_{k}\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}\otimes\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}$

11 for  _$n\leftarrow 1$ to $N$ _ do

$(u^{1,n},u^{2,n})^{T}\sim\mathcal{N}(\hat{\mu},\hat{\Sigma})$ // sample from
posterior

12
$G^{n,j}\leftarrow\sum_{k=1}^{J}u^{1,n}_{k}\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}$

13
$H^{n,j}\leftarrow\sum_{k=1}^{J}u^{2,n}_{k}\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}\otimes\frac{x^{k}-x^{j}}{\|x^{k}-x^{j}\|}$

14

Algorithm 2 Bayesian Ensemble-based Gradient Inference

######  Remark 2.6.

The numerical complexity of Algorithms 1 and 2 is dominated by the cost of
solving a linear equation of size $J\times J$, since we assume pointwise
evaluation of $V(x^{i})$ to be available initially.

###  2.2 Examples

We illustrate the proposed approximation of gradient and Hessian with the help
of two well-known benchmarks for optimization.

#### Rastrigin function in 1d

We set $V:\mathbb{R}\to\mathbb{R}$ with $V(x)=x^{2}+3\cdot(1-\cos(2\pi x))$,
which is (a slightly differently scaled variant of) the Rastrigin function.
Figure 1 demonstrates how we can locally approximate a function from given
pointwise evaluations. Note in particular how $\xi$ can be used as a parameter
for the degree of (non-)locality of the gradient approximation: $\xi=0$
corresponds to local Taylor approximation while the case $\xi\gg 1$ approaches
a least-squares quadratic fit through all data points (albeit always going
through the current reference point $(x^{j},V(x^{j}))$). The figure on the
bottom right also shows local quadratic approximation samples obtained by
interpreting (2.2) as a Bayesian inverse problem (after defining a prior on
the coefficients $u^{j}$), taking i.i.d. samples $\hat{u}^{j}$ from the
posterior. All local approximate functions are built as follows: Starting from
a coefficient vector ${\bar{u}}^{j}$ (i.e. either least squares solution to
(2.2) obtained with Algorithm 1 or posterior samples from Algorithm 2) we
construct a quadratic approximative map via

$x\mapsto V(x^{j})+\left\langle
G^{j}_{\xi}(\\{x^{i}\\}_{i=1}^{J}),~{}x-x^{j}\right\rangle+\frac{1}{2}\left\langle
x-x^{j},~{}H^{j}_{\xi}(\\{x^{i}\\}_{i=1}^{J})\cdot(x-x^{j})\right\rangle,$

and plot this quadratic function alongside the graph of $V$.

(a) Least squares solutions, $\xi=0$

(b) Least squares solutions, $\xi=1$

(c) Least squares solutions, $\xi=1000$

(d) $\xi=0$ and quadratic approximations obtained from samples from posterior
distribution on coefficients $u^{j}$.

Figure 1: Approximation of (a variant of) the Rastrigin function $x\mapsto
x^{2}+3(1-\cos(2\pi x))$ via finite evaluations with various choices of $\xi$.
Code: 1d_gradinf.py

#### Himmelblau function in 2d

As $V$ we choose the Himmelblau function
$V:(x,y)\mapsto(x^{2}+y-11)^{2}+(x+y^{2}-7)^{2}$ on $\mathbb{R}^{2}$ and
proceed along the same lines as for the Rastrigin function. Figure 2 shows the
influence of $J$ and $\xi$. The column on the right-hand side shows the level
sets of

$x\mapsto
V(x^{j})+\left(G^{j}(\\{x^{i}\\}_{i=1}^{J})\right)^{T}(x-x^{j})+\frac{1}{2}(x-x^{j})^{T}\cdot
H^{j}(\\{x^{i}\\}_{i=1}^{J})\cdot(x-x^{j})$ (2.5)

for an arbitrarily chosen “reference particle” $x^{j}$.

The left column of Figure 2 shows (for an ensemble $(x^{j},V(x^{j}))_{j}$) a
comparison between true gradients $\nabla V(x^{j})$ and the least squares
approximations obtained from Algorithm 1. Note that the arrows are scaled to
the same magnitude, so only the angle between ground truth and inferred
gradient is a relevant measure of approximation quality. Note that taking a
larger ensemble improves the local gradient approximation, as to be expected.
The right column shows level sets of the quadratic function in Equation 2.5.
For $\xi=0$, the quadratic approximation to $V$ obtained from (inferred)
gradient and Hessian works well only locally, while the global approximation
$\xi=1000$ shows deficits when judging actual gradient approximation (as can
be seen by comparing the arrows in the figures belonging to $J=25,\xi=0$ with
$J=25,\xi=1000$). This is a similar phenomenon as shown in Figure 1. Finally,
Figure 3 shows an illustration of Algorithm 2. It shows how the true gradient
$\nabla V(x^{i})$ in a specific ensemble member $x^{i}$ is approximated. In
the smaller ensemble of size $J=5$, the posterior gradient samples
predominantly point along a direction different to the true gradient, but the
latter can be seen to be in the support of this measure. By increasing the
ensemble size to $J=25$ and thereby providing more information, it can be seen
that the support of the posterior measure aligns with and contracts on the
true gradient.

(a) $J=5$, $\xi=0$

(b) quadratic approximation centered at $x^{j}$ (marked with square).

(c) $J=25$, $\xi=0$

(d) quadratic approximation centered at $x^{j}$ (marked with square).

(e) $J=25$, $\xi=1000$

(f) quadratic approximation centered at $x^{j}$ (marked with square).

Figure 2: Approximation of Himmelblau function by pointwise evaluation. Left
column shows approximated gradients $G^{j}(\\{x^{i}\\},V)$ in comparison to
actual gradients (vectors are rescaled to unit length). Right columns show
quadratic function obtained by extending local quadratic inferred
approximation to whole domain. Code: 2d_gradinf_J5_xi0.py,
2d_gradinf_J25_xi0.py, 2d_gradinf_J25_xi1000.py

(a) $J=5$

(b) $J=25$

Figure 3: Posterior gradient samples after choosing a prior on the
coefficients $u^{j}$: Radial histogram of gradient (direction) samples,
evaluated at a specific ensemble member, with true gradient (blue arrow) for
reference. Clearly, increasing the number of ensemble members $J$ contracts
the Bayesian posterior around the ground truth gradient vector. Note that
magnitude of any gradient vectors is disregarded in this illustration. Code:
2d_bayes_gradinf.py

We now describe a series of examples of optimization and sampling algorithms
which can be augmented with the gradient approximation presented above.

##  3 Consensus-based optimization augmented by EGI

For optimization problems with unknown structure it is advantageous to have
gradient-free optimization algorithms that allow for treating the objective as
black-box. Therefore the original formulation of CBO does not use any gradient
information. With the gradient inference idea discussed above, we are able to
keep the gradient-free property of CBO and still augment the algorithm with
higher order differential information.

###  3.1 EGI-CBO

Let $V:\mathbb{R}^{d}\to[0,\infty)$ be a (possibly non-convex) function,
$\beta>0$ an inverse heat parameter. We propose the following gradient-EGI-CBO
method describing the dynamics of the particles via

$\displaystyle\,\mathrm{d}x^{i}=-\kappa
G^{0}_{\xi}(\\{\mathsf{m}_{\beta}(\rho),x^{1},\ldots,x^{J}\\})\,\mathrm{d}t-\lambda(x^{i}-\mathsf{m}_{\beta})\,\mathrm{d}t+\sigma|x^{i}-\mathsf{m}_{\beta}|\,\mathrm{d}W_{t}^{i}$
(3.1)

with initial conditions $x_{0}^{i}\sim\mathcal{P}_{2}(\mathbb{R}^{d})$ drawn
independently. Note that the ensemble is augmented by the weighted mean and we
take the least square gradient and Hessian approximation centered at the
weighted mean and use this gradient approximation for all particles. The full
algorithm for this method – EGI-CBO – can be found in Algorithm 3.

The difference to CBO is the additional (approximated) gradient term. The
motivation for incorporating this is threefold: (1) Assuming that the spatial
resolution of the ensemble is coarser than the distance of individual minima
and the ensemble is still spread out, the gradient term can help to jump over
local minima. (2) The gradient term converges to the true gradient projected
onto the affine subspace spanned by the ensemble. This may facilitate the
convergence to the true global minimizer. Indeed, by the Laplace principle and
the quantitative non-asymptotic Laplace principle Fornasier et al. (2021b) the
approximation quality of CBO strongly depends on the temperature parameter.
Additional gradient information is expected to improve the approximation
quality, and drive the weighted mean towards the true minimizer. (3) In the
later phase of the dynamics, the additional gradient term helps with
accelerated collapse around the actual position of the (possibly local)
minimum. This alleviates an issue of CBO where the ensemble does indeed
collapse in a vicinity of a (possibly local) minimum $x^{\star}$, but with
consensus not necessarily converging to the actual position of $x^{\star}$.
Both features, “convergence to better minima” and “better convergence to
minima”, can be seen from the following numerical experiments.

###  3.2 Numerical examples: EGI-CBO

For simplicity, we use the least squares version (Algorithm 1) for the
gradient approximation for the following results. Note that a sampling
approach via Algorithm 2 is feasible as well and might be preferred in Machine
Learning applications as the behaviour resembles stochastic gradient descent
methods.

Data:
$N\in\mathbb{N},\\{x_{0}^{j}\\}_{j=1}^{J},\alpha,\lambda,\sigma,\kappa,\xi,\tau\geq
0$, noise$~{}\in\\{\text{norm-proportional},\text{component-wise}\\}$

Result: minimizer $\mathsf{m}_{\beta}$ of $I$

1 for  _$n\leftarrow 0$ to $N-1$ _ do

$m_{n}^{\alpha}\leftarrow\frac{\sum_{j=1}^{J}\exp(-\alpha I(x_{n}^{j}))\cdot
x_{n}^{j}}{\sum_{j=1}^{J}\exp(-\alpha I(x_{n}^{j}))}$ // Weighted mean -- use
logsumexp to avoid underflow

$m_{n}\leftarrow\frac{1}{J}\sum_{j=1}^{J}x_{n}^{j}$ // Unweighted mean

2 for  _$j\leftarrow 1$ to $J$ _ do

$g_{n}\leftarrow G^{0}(\\{m_{n}\\}\cup\\{x_{n}^{i}\\}_{i=1}^{J},I)$ //
Gradient approximation (least squares) at unweighted mean, via Algorithm 1

3 $W_{n}^{j}\sim N(0,I^{d})$

4 switch  _noise_ do

5 case  _norm-proportional_ do

6 $\sigma_{n}^{j}\leftarrow\sigma\|x_{n}^{j}-m_{n}^{\alpha}\|W_{n}^{j}$

7 case  _component-wise_ do

8 $\sigma_{n}^{j}\leftarrow\sigma(x_{n}^{j}-m_{n}^{\alpha})\odot W_{n}^{j}$

$x_{n+1}^{j}\leftarrow x_{n}^{j}-\tau\kappa
g_{n}-\tau\lambda(x_{n}^{j}-m_{n}^{\alpha})+\sqrt{\tau}\sigma^{j}_{n}$ //
Euler-Maruyama step

9

Algorithm 3 EGI-CBO

######  Remark 3.1.

If we choose $\lambda=0$ and $\sigma=0$ in Algorithm 3, then this corresponds
to a kind of ensemble-based gradient descent method. Unfortunately this seems
to perform quite badly in practice: When the ensemble moves to a position such
that its least squares approximation has a vanishing gradient (and this
happened frequently in our numerical experiments), then the iteration
essentially stops, without approaching the minimizer any further. For example,
imagine $V(x)=x^{2}$, and ensemble of size $J=2$ at positions $-1$ and $+1$.
Then the estimated gradient will be $0$, i.e. the ensemble will not move any
further (although this issue can be alleviated by extrapolating gradients via
second-order differential information, as done in Algorithm 7). The impact of
contraction and stochastic exploration means that this cannot happen easily
with EGI-CBO.

#### Rastrigin function in 2d

We aim to minimize the two-dimensional Rastrigin function
$V(x_{1},x_{2})=2\cdot 10+x_{1}^{2}-10\cos(2\pi x_{1})+y_{1}^{2}-10\cos(2\pi
y_{1})$. We start with an ensemble of $J=4$ particles, which is a moderately
large ensemble in two dimensions, and we also consider the case of $J=2$,
which yields an ensemble spanning an affine subspace of dimension 1 in each
iteration. We further set $\alpha=100$, $\lambda=1.5$, $\sigma=0.7$,
$\kappa=0.5$, $\xi=0$, $\tau=0.01$ and $N=1000$ iterations. Note that we
choose the initial ensemble uniform in the set $[-4,-1]^{2}$, which does not
contain the global minimum. As discussed in Kalise et al. (2022), this is a
much harder case where only few studies of the CBO method exist.

The results can be observed in Figure 4: Here we compare CBO and EGI-CBO with
ensemble sizes $J=4$ and $J=2$. We run 100 Monte Carlo simulations (with the
same initial ensemble but independently sampled noise in the iterations) for
each of these four settings and plot the position of the final iteration’s
weighted mean as a black dot in the function domain. We also superimpose a 2d
histogram to give a sense of the distribution of final weighted means of the
MC simulations. We can make out the following features:

  * •

Convergence towards local minima: It can be observed that the CBO dynamics
tends to terminate in a local minimum (and never in the global minimum), but
on a position with nonvanishing slope of $f$. This is due to the fact that the
weighted mean acts as a barrier to the left of the ensemble: The particle
closest to the weighted mean experience little to none drift or diffusion
itself, with approaching particles being attracted to it. This is resolved by
the approximate gradient term in EGI-CBO.

  * •

Exploration of “better minima”. Due to the gradient term in EGI-CBO, the
ensemble is able to experience a substantial shift along the essential slope
of $f$, even surpassing some intermittent local minima. CBO on the other hand
can only exhibit deterministic contraction towards $\mathsf{m}_{\beta}$ and
stochastic exploration (which ignores the shape of $f$). Note that this effect
is stronger for “global” EGI-CBO with $\xi>0$, see below.

  * •

This holds to some extent also in the underdetermined case $J=2$, although
convergence towards local minima is only slightly accelerated: With only two
live particles, gradient information will never point directly towards local
minimum.

######  Remark 3.2.

In addition to the observations discussed above, the EGI-CBO turns out to
reliably find various local minima of the highly multi-modal function within
$100$ Monte Carlo runs. This opens a new field of applications for CBO and
might be worth to be investigated in future work.

(a) Initial ensemble, $J=4$

(b) Initial ensemble, $J=2$

(c) 100 runs of CBO ($J=4$).

(d) 100 runs of CBO ($J=2$)

(e) 100 runs of EGI-CBO ($J=4$).

(f) 100 runs of EGI-CBO ($J=2$).

Figure 4:  Code: testMC_cbo_2d_rastrigin.py, testMC_cbo_2d_rastrigin_J2.py,
testMC_aug_cbo_2d_rastrigin.py,testMC_aug_cbo_2d_rastrigin_J2.py

#### Himmelblau function in 2d

In contrast to the Rastrigin function with steep valleys the Himmelblau
function is rather flat in the neighborhood of its minima. We also compare the
behaviour of CBO and EGI-CBO in this setting. The Himmelblau function exhibits
four global minima with value $0$. Figure 5 shows a comparison of a single run
of CBO with $64$ particles and a single run of EGI-CBO with only $3$
particles. We remark that all our experiments show similar results, so we show
only one specific instance. This allows us to better represent the dynamics.
It can be observed that EGI-CBO actually converges to one of the global
minima, with an exponential rate, while the ensemble of CBO collapses to a
point close to a minimum, but not towards the minimum itself. The function
value of $\mathsf{m}_{\beta}$ decreases in discrete jumps. This corresponds to
ensemble members crossing “better level sets” which influences
$\mathsf{m}_{\beta}.$ EGI-CBO on the other hand strongly profits from a
continuous drift towards a minimum.

(a) Initial ensemble $J=64$ (black dots) with initial weighted mean (star)

(b) One representative run of CBO ($J=64$).

(c) Initial ensemble $J=3$ (black dots) with initial weighted mean (star)

(d) One representative run of EGI-CBO ($J=3$).

Figure 5: Application of CBO and EGI-CBO to the Himmelblau test function.
Code: test_cbo_2d_himmelblau.py, test_aug_cbo_2d_himmelblau.py

#### High-dimensional smooth example

To further highlight the advantages of the gradient augmentation, we employ a
high-dimensional smooth convex example,
$V:x\mapsto\frac{1}{2}\|x-(1,\ldots,1)\|^{2}$ on $\mathbb{R}^{d}$, where
$d=10$. Clearly, $x=(1,\ldots,1)$ is the global minimum of $V$. We sample the
initial ensemble uniformly on $[-4,-1]^{d}$ and emphasize that this domain
excludes the global minimum. For the parameters, we set $\alpha=100$,
$\lambda=1$, $\sigma=0.2$ and we use component-wise noise for dimension-
robustness. We start with the over-determined case $J=20$. EGI-CBO
($\kappa=4.0$) exhibits exponential convergence, which is not that surprising
given that it essentially performs gradient descent, with its ensemble size
having enough descriptive power to span the full space. Vanilla CBO does not
converge to the minimum due to contraction of the ensemble on a point
different from the minimum. The fact that the initial ensemble is far away
from the minimum is an additional adversarial factor. Further experiments show
that this is not resolved by increasing the noise term as the stochastic
dynamics becomes unstable for a noise level larger than a certain threshold,
even for component-wise noise.

######  Remark 3.3.

We want to emphasize that this is not in contradiction with the theory for CBO
Carrillo et al. (2018); Fornasier et al. (2021b) as the proofs consider the
mean-field setting where the diffusion instantaneouly extends the support of
the distribution to the whole domain, and many proofs assume additionally that
the unique global minimizer is contained in the support of the initial
distribution.

The case $J=5$ is more interesting: While CBO performs badly (as to be
expected, see Remark 3.3), EGI-CBO is subjected to enough noise that
accumulated gradient information brings this underdetermined ensemble quite
close to the global minimum, although the ensemble has size lower than the
spatial dimension.

(a) CBO, $J=5$

(b) EGI-CBO, $J=5$

(c) CBO, $J=20$

(d) EGI-CBO, $J=20$

Figure 6: Performance of CBO and EGI-CBO on the function
$V:x\mapsto\|x-(1,\ldots,1)\|^{2}$ in $\mathbb{R}^{10}$. Diagrams show
evaluation of $V$ on the weighted mean over the course of the whole iteration.
Note that CBO is plotted over log-iterations in contrast to EGI-CBO, which is
plotted over linear iterations. CBO levels off close to the initialization,
with exponentially growing plateau lengths, while EGI-CBO converges
exponentially fast to the global minimum (more quickly for larger ensembles,
with jumps for smaller ensembles). Code: test_cbo_ndnorm.py

In the previous numerical examples the additional smoothing error term in
(2.1) is switched off ($\xi=0$). As discussed in Figure 1 this corresponds to
local gradient approximations. In next section we study the influence of
nonlocal gradient approximations.

###  3.3 Global versus Local EGI-CBO

Figure 1 showing global EGI for the one-dimensional Rastrigin function
suggests that “global” EGI-CBO (i.e. $\xi>0$) can lead to improved exploration
behaviour. In the following we investigate the influence of $\xi$ in more
detail. We therefore employ EGI-CBO as described in algorithm 3 with $\xi>0$.
This means that line 3 is modified to $g_{n}\leftarrow
G_{\xi}^{0}(\\{m_{n}\\}\cup\\{x_{n}^{i}\\}_{i=1}^{J})$ (i.e. with explicit
dependence on $\xi$), where we now compare the settings $\xi=0$ and $\xi>0$
for varying values. We set $\lambda=1.0$, $\sigma=0.5$ and in the augmented
case we have additionally $\kappa=2.0$ as well as $\xi=100$. Parameters of the
time discretization are $T=10$, $\tau=0.01$. In Figure 7 we note that
augmentation drives the ensemble towards the true minimizer. In contrast, the
reference solution of CBO gets stuck in higher level set regions.

(a) Initial ensemble ($J=3$) chosen for all simulations. Star marks weighted
mean.

(b) Histogram of final weighted mean of 100 runs of CBO. All results are equal
to the initial weighted mean of the ensemble.

(c) Histogram of final weighted mean of 100 runs of EGI-CBO ($\xi=0$).

(d) Histogram of final weighted mean of 100 runs of global EGI-CBO ($\xi>0$).

Figure 7: Demonstration of benefits of global EGI-CBO for multimodal
optimization. Code: testMC_cbo_globalapprox_1d_rastrigin.py

###  3.4 Discussion of EGI-CBO

To conclude the section on EGI-CBO we discuss some limitations of the method.

#### Subexponential convergence to higher-order minima.

For further investigation of the speed of convergence, we choose
$V(x)=\|x\|^{4}$ in $d=50$ and set $J=10$, $\kappa=2.5$, $\sigma=0.2$,
$\lambda=2.5$. In contrast to $\|x-(1,\ldots,1)\|^{2}$, where EGI-CBO
performed very nicely, we note that the gradient of $V$ vanishes rapidly in
the neighborhood of the global minimum vanishes. Figure 8 indicates that flat
regions around the minimizer slightly diminish the advantage of EGI-CBO over
vanilla CBO. Indeed, EGI-CBO exhibits only subexponential convergence.
Nevertheless, it still outperforms CBO which does not converge to the minimum
at all. (The latter fact is not demonstrated here, but experiments show the
same behaviour as for the test function $\|x-(1,\ldots,1)\|^{2}$ above).

Figure 8: Convergence toward a higher-order minimum.

#### Very high-dimensional multimodal functions.

In the literature Pinnau et al. (2017); Carrillo et al. (2021) very good
performance of CBO for the high-dimensional ($d=20$) Rastrigin function is
reported. In order to achieve this performance, it is necessary to either set
$J\gg d$ Pinnau et al. (2017) or $\sigma\gg\lambda$ Carrillo et al. (2021). In
the latter case the diffusion dominates the drift towards the weighted mean,
thus leaving the contractive domain of CBO which was established in the theory
(Carrillo et al., 2021, Theorem 3.2). In particular, CBO exhibits a random
exploration with a slight bias towards the weighted mean, and contraction of
the ensemble does not happen. Employing additional gradient information in
this non-contractive domain of CBO, leads to numerical instabilities which do
not improve the performance of vanilla CBO.

#### Smooth unimodal objective functions.

If $V$ is a smooth unimodal function with convex (or “nearly convex”)
structure – for example the typical inverse elliptical problem for strongly
informative data – CBO should not be used, as other methods typically have
much better performance. Similarly, if true gradient is cheap to obtain, CBO
is probably not the method of choice as it requires many function evaluations.
Nethertheless, for real-world applications often no structure is known
a-priori and only functions evaluations are available. Then (augmented) CBO is
a competitive alternative to other population based methods, which is
furthermore backed-up by theoretical convergence proofs.

##  4 Gradient-based sampling using EGI

In contrast to our approach with CBO (a completely gradient-agnostic method
which we augmented via inexact gradients), we now describe how some gradient-
based sampling methods can be adapted to inexact ensemble-based gradients.

###  4.1 EGI-LS and EGI-MALA

We begin with the unconditioned case, i.e. $M=\mathrm{Id}$ and propose the
following gradient-free sampling method as substitute dynamics:

$\displaystyle\,\mathrm{d}x^{j}_{t}=-G^{j}(\\{x^{i}\\}_{i=1}^{J})\,\mathrm{d}t+\sqrt{2}\,\mathrm{d}W^{j}_{t}$
(4.1)

with initial condition $x_{0}^{i}\sim\mathcal{P}_{2}(\mathbb{R}^{d})$
independently drawn for $i=1,\dots,J.$ A pseudocode for the implementation of
the dynamics is given in Algorithm 4.

In the same spirit, we provide results for the other ensemble-based sampling
algorithms presented in the introduction (see Section 1.1). Note that we set
$\xi=0$ by default. We expect that global gradient approximations lead to
similar behaviour as in the optimization case. For some sampling applications
it might make sense to allow for “global gradient approximations” (i.e.
$\xi>0$). For the sake of legibility we leave a detailed study to future
investigations.

An Metropolis adjustment, the so-called MALA sampler, was proposed to improve
the sampling accuracy of discretized Langevin dynamics in Roberts and Tweedie
(1996). In our setting this yields a straight-foward generalization of MALA,
the EGI-augmented Metropolis-adjusted Langevin algorithm (EGI-MALA), see
Algorithm 5.

Data: $N\in\mathbb{N},\\{x_{0}^{j}\\}_{j=1}^{J}$, step size $h$

Result: samples $\\{x_{n}^{j}\\}_{j=1,n=1}^{J,N}$ from
$\,\mathrm{d}\mu/\,\mathrm{d}x\propto\exp(-V(x))$

1 for  _$n\leftarrow 0$ to $N-1$ _ do

2 for  _$j\leftarrow 1$ to $J$ _ do

$g_{n}^{j}=G^{j}(\\{x_{n}^{i}\\}_{i=1}^{J})$// approximate gradient in
proposals, via Algorithm 1

3 $W_{n}^{j}\sim N(0,1)$

4 $x_{n+1}^{j}\leftarrow x_{n}^{j}-hg_{n}^{j}+\sqrt{2h}W^{j}_{n}$

5

6

Algorithm 4 EGI-augmented Langevin Sampler (EGI-LS)

Data: $N\in\mathbb{N},\\{x_{0}^{j}\\}_{j=1}^{J}$, $\tau>0$

Result: samples $\\{x_{n}^{j}\\}_{j=1,n=1}^{J,N}$ from
$\,\mathrm{d}\mu/\,\mathrm{d}x\propto\exp(-V(x))$

1 for  _$n\leftarrow 0$ to $N-1$ _ do

2 for  _$j\leftarrow 1$ to $J$ _ do

3 if  _$n\geq 1$_ then

$g_{n}^{j}=G^{j}(\\{x_{n}^{i}\\}_{i=1}^{J}\cup\\{\mu_{n-1}^{i}\\}_{i=1}^{J})$
// approximate gradient, via Algorithm 1

4

5 else

$g_{n}^{j}=G^{j}(\\{x_{n}^{i}\\}_{i=1}^{J})$ // approximate gradient, via
Algorithm 1

6

7 $W_{n}^{j}\sim N(0,I^{d})$

$\mathrm{prop}_{n}^{j}\leftarrow x_{n}^{j}-\tau
g_{n}^{j}+\sqrt{2\tau}W^{j}_{n}$// Langevin-type proposal

$\gamma_{n}^{j}=G^{j}(\\{\mathrm{prop}_{n}^{i}\\}_{i=1}^{J})$// approximate
gradient in proposals, via Algorithm 1

8

9 for  _$j\leftarrow 1$ to $J$ _ do

$q_{\mathrm{fwd}}^{j}\leftarrow\exp\left(-\frac{1}{4\tau}\|\mathrm{prop}_{n}^{j}-(x_{n}^{j}-\tau
g_{n}^{j})\|^{2}\right)$// Metropolis adjustment

10

11
$q_{\mathrm{bwd}}^{j}\leftarrow\exp\left(-\frac{1}{4\tau}\|x_{n}^{j}-(\mathrm{prop}_{n}^{j}-\tau\gamma_{n}^{j})\|^{2}\right)$

12 $\alpha_{j}\leftarrow\frac{\exp(-V(\mathrm{prop}_{n}^{j}))\cdot
q_{\mathrm{bwd}}^{j}}{\exp(-V(x_{n}^{j}))\cdot q_{\mathrm{fwd}}^{j}}$

13 $\xi\sim\mathrm{Unif}[0,1]$

14 if  _$\xi\leq\alpha_{j}$_ then

15 $x_{n+1}^{j}\leftarrow\mathrm{prop}_{n}^{j}$

$\mu_{n}^{j}\leftarrow x_{n}^{j}$ // keep previous iterate as memory

16

17 else

18 $x_{n+1}^{j}\leftarrow x_{n}^{j}$

$\mu_{n}^{j}\leftarrow\mathrm{prop}_{n}^{j}$ // Keep (rejected) proposal as
memory

19

Algorithm 5 (EGI-MALA)

###  4.2 EGI-ALDI and EGI-EKS

The extraordinary feature of gradient-free ALDI is that it implicitly
constructs and uses gradient information of $V$ without the need of solving a
linear system as demonstrated in 1.1. Unfortunately, this approximation fails
if $A$ is nonlinear.

For this reason, we also consider a new version of this methodology by
replacing gradient-free ALDI’s implicit gradient (biased for nonlinear forward
maps) by our inexact gradient approximation via EGI, leading to the following
sampling algorithm, EGI-ALDI (gradient-free ALDI with estimated gradients):

Data: $N\in\mathbb{N},\\{x_{0}^{j}\\}_{j=1}^{J}$, $\tau>0$

Result: samples $\\{x_{n}^{j}\\}_{j=1,n=1}^{J,N}$ from
$\,\mathrm{d}\mu/\,\mathrm{d}x\propto\exp(-V(x))$

1 for  _$n\leftarrow 0$ to $N-1$ _ do

2 $\bar{x}_{n}\leftarrow\frac{1}{J}\sum_{i=1}^{J}x_{n}^{i}$

3
$C_{n}^{1/2}\leftarrow\frac{1}{\sqrt{J}}\left(x_{n}^{1}-\bar{x}_{n},\ldots,x_{n}^{J}-\bar{x}_{n}\right)\in\mathbb{R}^{d\times
J}$

4
$C_{n}\leftarrow\frac{1}{J}\sum_{i=1}^{J}(x_{n}^{i}-\bar{x})\otimes(x_{n}^{i}-\bar{x})$

5 for  _$j\leftarrow 1$ to $J$ _ do

$g_{n}^{j}\leftarrow G^{j}(\\{x_{n}^{i}\\}_{i=1}^{J}\cup\\{\bar{x}_{n}\\})$ //
approximate gradient in each particle (not in mean), via Algorithm 1

6 $W^{j}_{n}\sim N(0,I^{J})$

$x_{n+1}^{j}\leftarrow x_{n}^{j}-\tau C_{n}\cdot
g_{n}^{j}+\tau\frac{d+1}{J}(x_{n}^{j}-\bar{x}_{n})+\sqrt{2\tau}C^{1/2}_{n}W^{j}_{n}$
// ALDI step

7

Algorithm 6 gradient-free ALDI with estimated gradients (EGI-ALDI)

###  4.3 Gradient extrapolation

All methods presented so far employ gradient approximations computed in each
ensemble point, see e.g. the for-loop in algorithm 4–4, which means that each
iteration requires the solution of $J$ linear equations (by evaluating
$G^{j}(\\{x_{n}^{i}\\}_{i},I)$ for $j=\\{1,\ldots,J\\}$). It is possible to
strongly cut down on computational complexity at the cost of introducing an
additional source of error: Instead of approximating each gradient, we can
approximate the gradient and Hessian only in a suitably chosen reference point
$x^{\star}$ (e.g. the ensemble mean) and extrapolate:

$\nabla V(x_{n}^{j})\approx
G^{0}(\\{x^{\star}\\}\cup\\{x_{n}^{i}\\}_{i=1}^{J})+H^{0}(\\{x^{\star}\\}\cup\\{x_{n}^{i}\\}_{i=1}^{J})\cdot(x_{n}^{j}-x^{\star}),$
(4.2)

which is reminiscent of $\nabla V(y)=\nabla V(x)+HV(x)\cdot(y-x)+o(\|y-x\|)$.
Of course, this approximation is only valid if either the ensemble is
sufficiently strongly concentrated or the measure is sufficiently Gaussian-
like. To be more precise, we describe this variant for the EGI-ALDI algorithm:

Data: $N\in\mathbb{N},\\{x_{0}^{j}\\}_{j=1}^{J}$, $\tau>0$

Result: samples $\\{x_{n}^{j}\\}_{j=1,n=1}^{J,N}$ from
$\,\mathrm{d}\mu/\,\mathrm{d}x\propto\exp(-V(x))$

1 for  _$n\leftarrow 0$ to $N-1$ _ do

2 $\bar{x}_{n}\leftarrow\frac{1}{J}\sum_{i=1}^{J}x_{n}^{i}$

3
$C_{n}^{1/2}\leftarrow\frac{1}{\sqrt{J}}\left(x_{n}^{1}-\bar{x}_{n},\ldots,x_{n}^{J}-\bar{x}_{n}\right)\in\mathbb{R}^{d\times
J}$

4
$C_{n}\leftarrow\frac{1}{J}\sum_{i=1}^{J}(x_{n}^{i}-\bar{x})\otimes(x_{n}^{i}-\bar{x})$

$g_{n}\leftarrow G^{0}(\\{\bar{x}_{n}\\}\cup\\{x_{n}^{i}\\}_{i=1}^{J})$ //
approximate gradient only in mean, via Algorithm 1

$H_{n}\leftarrow H^{0}(\\{\bar{x}_{n}\\}\cup\\{x_{n}^{i}\\}_{i=1}^{J})$ //
Hessian approximation in mean, via Algorithm 1

5 for  _$j\leftarrow 1$ to $J$ _ do

6 $g_{n}^{j}\leftarrow g_{n}+H_{n}(x_{n}^{j}-\bar{x}_{n})$

$W^{j}_{n}\sim N(0,I^{J})$// extrapolate estimated gradient

$x_{n+1}^{j}\leftarrow x_{n}^{j}-\tau C_{n}\cdot
g_{n}^{j}+\tau\frac{d+1}{J}(x_{n}^{j}-\bar{x}_{n})+\sqrt{2\tau}C^{1/2}_{n}W^{j}_{n}$//
ALDI step

7

Algorithm 7 gradient-free ALDI with estimated and extrapolated gradients (EGI-
ALDI-extra)

By taking EGI-ALDI-extra and using a sampled gradient approximation (Algorithm
2) instead of the least square solution (Algorithm 1) we obtain Algorithm 8,
although honestly mainly conceived for the sake of its acronym.

Data: $N\in\mathbb{N},\\{x_{0}^{j}\\}_{j=1}^{J}$, $\tau>0$

Result: samples $\\{x_{n}^{j}\\}_{j=1,n=1}^{J,N}$ from
$\,\mathrm{d}\mu/\,\mathrm{d}x\propto\exp(-V(x))$

1 for  _$n\leftarrow 0$ to $N-1$ _ do

2 $\bar{x}_{n}\leftarrow\frac{1}{J}\sum_{i=1}^{J}x_{n}^{i}$

3
$C_{n}^{1/2}\leftarrow\frac{1}{\sqrt{J}}\left(x_{n}^{1}-\bar{x}_{n},\ldots,x_{n}^{J}-\bar{x}_{n}\right)\in\mathbb{R}^{d\times
J}$

4
$C_{n}\leftarrow\frac{1}{J}\sum_{i=1}^{J}(x_{n}^{i}-\bar{x})\otimes(x_{n}^{i}-\bar{x})$

$g_{n}\leftarrow G^{0}(\\{\bar{x}_{n}\\}\cup\\{x_{n}^{i}\\}_{i=1}^{J})$ //
approximate gradient only in mean, via Algorithm 2

$H_{n}\leftarrow H^{0}(\\{\bar{x}_{n}\\}\cup\\{x_{n}^{i}\\}_{i=1}^{J})$ //
Hessian approximation in mean, via Algorithm 2

5 for  _$j\leftarrow 1$ to $J$ _ do

$g_{n}^{j}\leftarrow g_{n}+H_{n}(x_{n}^{j}-\bar{x}_{n})$// extrapolate
estimated gradient

6 $W^{j}_{n}\sim N(0,I^{J})$

$x_{n+1}^{j}\leftarrow x_{n}^{j}-\tau C_{n}\cdot
g_{n}^{j}+\tau\frac{d+1}{J}(x_{n}^{j}-\bar{x}_{n})+\sqrt{2\tau}C^{1/2}_{n}W^{j}_{n}$//
ALDI step

7

Algorithm 8 gradient-free ALDI with Bayesian Gradients (ALDI-BaG)

###  4.4 Numerical example: Sampling from a two-dimensional non-Gaussian
measure

We consider the inverse problem

$y=A(x)+\varepsilon$

with $A:\mathbb{R}^{2}\to\mathbb{R}$, $A(x)=(x_{2}-2)^{2}-(x_{1}-3.5)-1$,
$y=0$, and $\varepsilon\sim N(0,\tfrac{1}{2^{2}})$. We set a Gaussian prior
$\mu_{0}=N(0_{2},\tau^{2}\cdot I)$ with $\tau=2$ on the unknown parameter
$x\in\mathbb{R}^{2}$. We write $\Phi(x)=\frac{1}{2\sigma^{2}}|y-A(x)|^{2}$,
$I(x)=\Phi(x)+\tfrac{1}{2\tau^{2}}\|x\|^{2}$. This leads to a Bayesian
posterior on the parameter $\mu^{y}$ of form

$\frac{\,\mathrm{d}\mu^{y}}{\,\mathrm{d}\mu_{0}}(x)\propto\exp\left(-\Phi(x)\right)$

or

$\frac{\,\mathrm{d}\mu^{y}}{\,\mathrm{d}x}(x)\propto\exp\left(-I(x)\right).$

Figure 9 shows the result of applying various samplers. We compare Ensemble
Langevin (EGI-LS), Ensemble MALA (EGI-MALA), gradient-free ALDI (Garbuno-Inigo
et al. (2020b), as an efficient gradient-free ensemble-based sampler),
gradient-free ALDI augmented by approximated gradient information (EGI-ALDI),
gradient-free ALDI with estimated and extrapolated gradients (EGI-ALDI-extra),
and the CBS sampler of Carrillo et al. (2022a). We start with a very small
ensemble of size $J=2$ and $N=10000$. It can be observed that EGI-LS needs a
larger number of ensemble members to perform well. This is due to the fact
that $J=2$ does not yield a sufficiently good approximation to the gradient.
Gradient-free ALDI, EGI-ALDI and EGI-ALDI-extra are restricted to the one-
dimensional affine subspace spanned by the initial ensemble and also do not
sample correctly from the posterior. CBS cannot extract sufficient information
from its two-point ensemble: Both the subspace property as well as the fact
that (in contrast to ALDI) no implicit gradient information is obtained leads
to quick collapse to a point in the vicinity of the “better” ensemble member.
EGI-MALA yields good samples: The shortcomings of EGI-LS are reigned in by the
Metropolis adjustment (atrociously bad proposals due to faulty gradient
approximations are rejected). For ensemble size $J=20$ (and $N=2000$) we see
that all methods drastically improve: EGI-LS can improve its fit with the
measure due to better gradient approximation (but still shows some bias), and
ALDI is no longer constrained to a lower-dimensional affine subspace.
Nevertheless, the forward mapping is sufficiently nonlinear that gradient-free
ALDI incorporates substantial bias which strongly impacts its performance. CBS
suffers in a similar way from its restricted applicability to non-Gaussian
measures. EGI-ALDI and EGI-ALDI-extra profit from the gradient approximation
but do not match EGI-MALA, which provides a near-perfect representation of the
measure. Of course, this comes at the cost of needing to solve a larger number
of linear equations as the following list (showing the amount of additional
linear equations to solve per iteration) illustrates:

  * •

EGI-LS: $J$ linear equations (gradient approximation in each particle)

  * •

EGI-MALA: $2\cdot J$ linear equations (gradient approximation in each particle
and each proposal)

  * •

CBS: no additional linear equations

  * •

ALDI: no additional linear equations

  * •

EGI-ALDI: $J$ linear equations (gradient approximation in each particle)

  * •

EGI-ALDI-extra: $1$ linear equation (gradient approximation only in ensemble
mean)

######  Remark 4.1.

It is striking that EGI-ALDI-extra shows such a strong improvement over ALDI
by use of just one additional linear equation solve per iteration, with EGI-
ALDI not being much better (but much more expensive). EGI-MALA has the best
performance but requires the highest number of linear equations to solve.

(a) 2d histograms of samplers. $J=2$, $N=10000$

(b) 2d histograms of samplers. $J=20,N=2000$

(c) Marginal histogram of all six methods in comparison, where $J=2$,
$N=10000$

(d) Marginal histogram of all six methods in comparison, where $J=20,N=2000$

Figure 9: Comparison of EGI-LS, EGI-MALA, gradient-free ALDI and EGI-ALDI.
$J=2$, $N=10000$ (left column) and $J=20$, $N=4000$ (right column). Bottom row
(figures 9(c), 9(d)) show histograms of $x$-component of samples. The orange
line marks the true marginal distribution. Code: sampling_2d_J2.py,
sampling_2d_J20.py

###  4.5 Discussion of results

  * •

If the system is “almost linear”, gradient-free ALDI and CBS is much more
efficient and should be preferred over EGI-LS or EGI-MALA, as there is no need
to solve additional linear systems.

  * •

On the other hand, if the ensemble size is limited (for example if evaluation
of $V$ is expensive, e.g. if it involves the numerical solution of expensive
computational models), then gradient approximation is a cost-effective add-on
to completely gradient-agnostic methods like CBS and Random Walk Metropolis
Hastings.

  * •

If the measure is multimodal or has a non-negligibly curved shape and gradient
information is unavailable, EGI-LS and EGI-MALA can accurately approximate
measures.

  * •

If the amount of ensemble members is chosen relatively small ($J\approx d$ or
less), then EGI-MALA performs much better than EGI-LS due to the Metropolis-
adjustment involved.

  * •

The extrapolation idea can be used to cut down on the amount of linear
equations that needs to be solved, from $\mathcal{O}(J)$ to $\mathcal{O}(1)$.

##  5 Conclusion and outlook

We have described a way to turn implicit differential information from
pointwise ensemble evaluation into an estimator for first (and higher order)
derivatives via Ensemble-based Gradient Inference (EGI). We have presented a
novel way to augment Consensus-based optimization by adding an additional
drift term proportional to this approximated gradient term which can be used
to find better minima and accelerates local convergence. We have also
demonstrated that sampling algorithms which usually work with exact gradients
can be used with inexact gradients via EGI with similar performance, but
without the need for explicit gradient evaluation.

There is a lot of perspective for future work: More theoretical analysis could
shed some light on the right balancing between the three terms in EGI-CBO
(gradient, consensus-building, and exploratory diffusion). At this point it is
still unclear under which conditions EGI-CBO can be proven to be stable and/or
what its limit points are. Similarly, the heuristical trade-off between local
($\xi=0$) and global ($\xi>0$) EGI-CBO would need to be studied in more depth
and this pertains to the sampling methods presented as well. In optimization
settings where there are several similar optima one might be interested in
finding all of them. A localization of CBO could help solve this issue and is
subject to future work. Also for sampling methods, the analysis of convergence
properties will help to guide the development of EGI variants, in particular
the analysis of computational costs of EGI variants vs. accuracy improvement.

## Acknowledgements

PW acknowledges support from MATH+ project EF1-19: Machine Learning Enhanced
Filtering Methods for Inverse Problems.

## References

  * Back (1996) T. Back.  _Evolutionary algorithms in theory and practice: evolution strategies, evolutionary programming, genetic algorithms_.  Oxford University Press, 1996.
  * Blömker et al. (2019) D. Blömker, C. Schillings, P. Wacker, and S. Weissmann.  Well posedness and convergence analysis of the ensemble Kalman inversion.  _Inverse Problems_ , 35(8):085007, 2019.
  * Borghi et al. (2022) G. Borghi, M. Herty, and L. Pareschi.  An adaptive consensus based method for multi-objective optimization with uniform Pareto front approximation.  _arXiv:2208.01362_ , 2022.
  * Carrillo et al. (2022a) J. Carrillo, F. Hoffmann, A. Stuart, and U. Vaes.  Consensus-based sampling.  _Studies in Applied Mathematics_ , 148(3):1069–1140, 2022a.
  * Carrillo et al. (2018) J. A. Carrillo, Y.-P. Choi, C. Totzeck, and O. Tse.  An analytical framework for consensus-based global optimization method.  _Mathematical Models and Methods in Applied Sciences_ , 28(06):1037–1066, 2018.
  * Carrillo et al. (2021) J. A. Carrillo, S. Jin, L. Li, and Y. Zhu.  A consensus-based global optimization method for high dimensional machine learning problems.  _ESAIM: Control, Optimisation and Calculus of Variations_ , 27:S5, 2021.
  * Carrillo et al. (2022b) J. A. Carrillo, F. Hoffmann, A. M. Stuart, and U. Vaes.  Consensus-based sampling.  _Studies in Applied Mathematics_ , 148(3):1069–1140, 2022b.
  * Fornasier et al. (2020) M. Fornasier, H. Huang, L. Pareschi, and P. Sünnen.  Consensus-based optimization on hypersurfaces: Well-posedness and mean-field limit.  _Math. Mod. Meth. Appl. Sci._ , 30(14):2725–2751, 2020.
  * Fornasier et al. (2021a) M. Fornasier, H. Huang, L. Pareschi, and P. Sünnen.  Consensus-based optimization on the sphere: Convergence to global minimizers and machine learning.  _Journal of Machine Learning Research_ , 22:1–55, 2021a.
  * Fornasier et al. (2021b) M. Fornasier, T. Klock, and K. Riedl.  Consensus-based optimization methods converge globally.  _arXiv:2103.15130_ , 2021b.
  * Fornasier et al. (2022) M. Fornasier, T. Klock, and K. Riedl.  Convergence of anisotropic consensus-based optimization in mean-field law.  _Applications of Evolutionary Computation_ , 13224:738–754, 2022.  doi: 10.1007/978-3-031-02462-7_46.
  * Garbuno-Inigo et al. (2020a) A. Garbuno-Inigo, F. Hoffmann, W. Li, and A. M. Stuart.  Interacting langevin diffusions: Gradient structure and ensemble kalman sampler.  _SIAM Journal on Applied Dynamical Systems_ , 19(1):412–441, 2020a.
  * Garbuno-Inigo et al. (2020b) A. Garbuno-Inigo, N. Nüsken, and S. Reich.  Affine invariant interacting Langevin dynamics for Bayesian inference.  _SIAM Journal on Applied Dynamical Systems_ , 19(3):1633–1658, 2020b.
  * Grassi and Pareschi (2021) S. Grassi and L. Pareschi.  From particle swarm optimization to consensus based optimization: Stochastic modeling and mean-field limit.  _Math. Mod. Meth. Appl. Sci._ , 31(08):1625–1657, 2021.
  * Ha et al. (2022) S.-Y. Ha, M. Kang, D. Kim, J. Kim, and I. Yang.  Stochastic consensus dynamics for nonconvex optimization on the stiefel manifold: Mean-field limit and convergence.  _Math. Mod. Meth. Appl. Sci._ , 32(03):533–617, 2022.
  * Huang (2021) H. Huang.  A note on the mean-field limit for the particle swarm optimization.  _Applied Mathematics Letters_ , 117:107133, 2021.
  * Huang and Qiu (2022) H. Huang and J. Qiu.  On the mean-field limit for the consensus-based optimization.  _Mathematical Methods in the Applied Sciences_ , 45(12):7814–7831, 2022.
  * Jordan et al. (1998) R. Jordan, D. Kinderlehrer, and F. Otto.  The variational formulation of the Fokker–Planck equation.  _SIAM Journal on Mathematical Analysis_ , 29(1):1–17, 1998.  URL https://doi.org/10.1137/S0036141096303359.
  * Kalise et al. (2022) D. Kalise, A. Sharma, and M. V. Tretyakov.  Consensus based optimization via jump-diffusion stochastic differential equations.  _arXiv:2205.04880_ , 2022.
  * Klamroth et al. (2022) K. Klamroth, M. Stiglmayr, and C. Totzeck.  Consensus-based optimization for multi-objective problems: a multi-swarm approach.  _in preparation_ , 2022.
  * Nüsken and Reich (2019) N. Nüsken and S. Reich.  Note on interacting Langevin diffusions: Gradient structure and ensemble Kalman sampler by Garbuno-Inigo, Hoffmann, Li and Stuart.  _arXiv preprint arXiv:1908.10890_ , 2019.
  * Pinnau et al. (2017) R. Pinnau, C. Totzeck, O. Tse, and S. Martin.  A consensus-based model for global optimization and its mean-field limit.  _Mathematical Models and Methods in Applied Sciences_ , 27(01):183–204, 2017.
  * Roberts and Tweedie (1996) G. O. Roberts and R. L. Tweedie.  Exponential convergence of Langevin distributions and their discrete approximations.  _Bernoulli_ , pages 341–363, 1996.
  * Schillings and Stuart (2017) C. Schillings and A. M. Stuart.  Analysis of the ensemble Kalman filter for inverse problems.  _SIAM Journal on Numerical Analysis_ , 55(3):1264–1290, 2017.
  * Simon (2013) D. Simon.  _Evolutionary optimization algorithms_.  John Wiley & Sons, 2013.
  * Totzeck (2022) C. Totzeck.  Trends in consensus-based optimization.  In N. Bellomo, J. A. Carrillo, and E. Tadmor, editors, _Active Particles, Volume 3: Advances in Theory, Models, and Applications_ , pages 201–226. Springer International Publishing, Cham, 2022.
  * Totzeck and Wolfram (2020) C. Totzeck and M.-T. Wolfram.  Consensus-based global optimization with personal best.  _Math. Biosci. Eng._ , 17(5):6026–6044, 2020\.  doi: 10.3934/mbe.2020320.
  * van Laarhoven and Aarts (1987) P. J. van Laarhoven and E. H. Aarts.  _Simulated Annealing: Theory and Applications_.  Springer, 1987.
  * Yang et al. (2019) J. Yang, G. O. Roberts, and J. S. Rosenthal.  Optimal scaling of Metropolis algorithms on general target distributions.  _arXiv: Computation_ , 2019.


2012 Vol. X No. XX, 000–000

11institutetext:  National Astronomical Observatories$/$Nanjing Institute of
Astronomical Optics $\&$ Technology, Chinese Academy of Sciences, Nanjing
210042,China<EMAIL_ADDRESS>
22institutetext: Key Laboratory of Astronomical Optics $\&$ Technology,
National Astronomical Observatories$/$Nanjing Institute of Astronomical Optics
$\&$ Technology, Chinese Academy of Sciences, Nanjing 210042,China
33institutetext: Physics $\&$ Astronomy Department, California State
University Northridge, 18111 Nordhoff Street, Northridge, California
91330-8268, USA
44institutetext: University of Chinese Academy of Sciences, Beijing 100049,
China

# A high-contrast coronagraph for earth-like exoplanet direct imaging: design
and test ∗ 00footnotetext: $*$ Supported by the National Natural Science
Foundation of China and Strategic Priority Research Program, CAS.

Cheng-Chao Liu 112244 De-Qing Ren 112233 Jian-Pei Dou 1122 Yong-Tian Zhu 1122
Xi Zhang 1122 Gang Zhao 1122 Zhen Wu 1122 Rui Chen 1122

###### Abstract

The high-contrast coronagraph for direct imaging earth-like exoplanet at the
visible needs a contrast of $10^{-10}$ at a small angular separation of
$4\lambda/D$ or less. Here we report our recent laboratory experiment that is
close to the limits. The test of the high-contrast imaging coronagraph is
based on our step-transmission apodized filter. To achieve the goal, we use a
liquid crystal array (LCA) as a phase corrector to create a dark hole based on
our dedicated focal dark algorithm. We have suppressed the diffracted and
speckle noise near the star point image to a level of 1.68$\times$$10^{-9}$ at
$4\lambda/D$, which can be immediately used for the direct imaging of Jupiter
like exoplanets. This demonstrates that high-contrast coronagraph telescope in
space has the potentiality to detect and characterize earth-like planets.

###### keywords:

instrumentation: coronagraph, liquid crystal array — techniques: apodization —
methods: laboratory

##  1 Introduction

Today, searching the universe for earth-like exoplanets is one of the hottest
topics in science. Up to now more than 1,000 exoplanets have been detected and
confirmed. Most of them were discovered through indirect ways such as the
transit method and radial velocity detection. These methods have not been
proven sensitive enough for the detection of earth-like planets. The Kepler
mission has discovered 2 earth-size planets: Kelper-62 in the habitable zone
and Kepler-78b by transit approach recently (Borucki et al. 2013;Pepe et al.
2013). However, direct imaging is needed to fully detect the signatures of
life or biomarkers $O_{2}$ and $O_{3}$ at visible wavelengths on the
exoplanets by using spectroscopic techniques. To realize this goal, the
unwanted diffraction light from host star and speckle noise from the optics
system must be efficiently suppressed using high-contrast imaging
coronagraphs. A contrast of $\sim 10^{-10}$ is required for imaging earth-like
planets in the visible wavelength, which is a great challenge for current
techniques. A critical issue is the speckle noise which is induced by any
optical imperfection, such as the wave-front error and surface roughness. The
static speckle noise, in generally degrades the contrast of a coronagraph to a
level on the order of $10^{-5}\sim 10^{-7}$ (Dou et al. 2010; Ren et al.
2010).

A common approach used to achieve this high-contrast by other research groups
is to correct the optical errors and generate a small local dark hole by using
one or two deformable mirrors (DMs), such as the groups led by Jet Propulsion
Laboratory (JPL) and Princeton (Trauger & Traub 2007; Kay et al. 2009).
However, this approach needs a deformable mirror (DM) with a large number of
actuators and complicated phase retrieval algorithm. Although the contrast can
reach the level of $10^{-10}$, the region of the dark holes created by this
two groups is not large enough to apply to space telescopes for exoplanet
detecting. In our group, we have successfully used the stochastic parallel
gradient descent (SPGD) algorithm to efficiently supress speckle noises
according to the point spread function (PSF) on the imaging focal plane, and
this algorithm is first used in our extreme adaptive optics (ExAO) for a
ground-based high-contrast system, which successfully demonstrated that a
small dark hole with an extra contrast gain of 65-time improvement is
achievable with a DM of 140 actuators (Ren et al. 2012). Recently, this
approach is used in our coronagraph, and the dark hole created by our approach
is much larger than that created by other groups, which we will discuss in
Section 3.

Ren and Zhu have developed a stepped transmission filter-based coronagraph
(Ren & Zhu 2007). Recently new designed filter can realize a contrast of
$10^{-6}$ at a small angular separation down to $4\lambda/D$ ($\lambda$ is the
wavelength and D is the diameter of telescope aperture) in our laboratory
test.

Based on the stepped filters, test results of the new designed coronagraph,
combined both the stepped transmission filters and our SPGD based dark-hole
correction, are presented in Section 2. In Section 3 we show the test optics
system and experimental results. Finally, we present the conclusions and
outline future developments.

##  2 TESTING RESULTS OF NEW STEPPED FILTERS AND THE DARK-HOLE CORRECTION
BASED ON SPGD ALGORITHM

###  2.1 New filters and the testing results

Ren $\&$ Zhu (2007) have firstly proposed a coronagraph based on stepped
transmission filters. This kind of coronagraph can be applied to an off-axis
space-based telescope which is aimed for a future space mission for direct
imaging of earth-like exoplanet. In our previous papers, we have demonstrated
this feasibility with a high-contrast imaging coronagraph (Dou et al. 2009).
Recently we have designed and manufactured five new filters with the goal for
better contrast over a large discovery area.

Here we briefly give the design and testing results in our laboratory. For
each filter, the diffraction is suppressed along one direction only (x or y),
depending on the orientation of the transmission attenuation in the filter. On
the coronagraph focal plane, the complex amplitude of the electric field is
the Fourier transform of the complex amplitude of the pupil. The theoretical
designed contrast for one filter should reach $10^{-6}$ at the angular
distance of equal or larger than $4\lambda/D$ from the PSF intensity peak. For
this design, two filters can also be combined together one after the other
with $90^{\circ}$ orientation. In this case diffraction along the two
perpendicular directions of x and y are suppressed, and a theoretical contrast
of $10^{-8}$ should be achieved at an angular distance of $4\lambda/D$ or
larger.

However, in practice the coronagraph can’t achieve the theoretical contrast.
Because the transmission of the filter is realized by uneven metallic coating,
which introduces a thickness variation of optical path, it results in an
unexpected phase. High contrast of this transmission filter based coronagraph
is achieved by the so-called apodized pupil. Due to the transmission apodized
pupil, the transmission throughput for each filter is 0.498. Figure 1 shows
the theoretical transmission patterns of the designed filters (top panel) and
the associated PSF images (bottom panels). The top panel shows images of one
filter (in the left-top panel), another filter with $90^{\circ}$ orientation
(in the central-top panel), and the combination of these two filters (in the
right-top panel). In figure 2, the theoretical stepped transmission of one
transmission filter and real photograph of the two orthogonal-orientation
filters are presented. PSF images of coronagraph based on the combination of
these two orthogonal-orientation filters under different exposure time of
0.15s, 15s and 1500s and associated contrast plot are shown in Figure 3. A
16-bit SBIG camera is used to take these images.

###  2.2 The dark-hole correction based on SPGD using LCA

Figure 1: Transmission patterns of the filters (top) and the associated PSFs
(bottom).

We have firstly proposed and demonstrated that SPGD algorithm can be applied
to the high-contrast coronagraph to reduce speckle noise induced by
imperfection optics (Dong et al. 2011). In previous paper, we used a
deformable mirror (DM) to provide an optimal phase to create a dark-hole in a
target region and the principle of SPGD algorithm has been discussed in
details somewhere (Ren et al. 2012). The DM used in his previous test has 140
effective actuators arranged in 12$\times$12 grid. The dark-hole has a small
area because of the small actuator numbers. Considering a DM (32$\times$32 or
64$\times$64 actuators) is too expensive to afford for our tests, in our
recent test, we use XY nematic series spatial light modulator (SLM) as a phase
corrector to replace a DM. The SLM is manufactured by BNS has 512$\times$512
pixels, with a 83.4% fill factor. In the test the pixels are divided into a
number of groups and each group is binned as an effective pixel consisting of
8$\times$8 pixels. The XY nematic series SLM is optimized to provide a full
wave ($2\pi$) of phase stroke upon reflection at the nominal design
wavelength.

In this test, we use SPGD optimization algorithm, which measures the wave-
front error according to the PSF on the focal plane, and the measured
information is used to directly control the phase of the LCA to create a dark-
hole.

Figure 2: The theoretical 51 step transmission filter and photograph of the
real two filters.

Figure 3: The top panels show the PSF images of the two filters coronagraph
under different exposure time. The bright vertical pattern in central and
right panels is due to the CCD image bloom. The bottom panel shows the
comparison between the test and designed PSF contrast along the diagonal
direction. Figure 4: The optical configuration of intensity calibration.

Figure 5: The intensity of starlight with the exposure of 0.5ms (left), the
intensity of planet light with the exposure of 85ms (right). Figure 6: The
whole optical system of dark-hole correction.

##  3 CONFIGURATIFig1ON OF DARK-HOLE CORRECTION OPTICS SYSTEM AND THE DARK
HOLE CREATION

###  3.1 The intensity calibration of starlight and planet light

Since a 16-bit CCD camera can only measure a contrast up to 1$/$65,536 (about
$10^{-4}$) and the contrast we will measure is out of this range, we need use
a neutral density (ND) filter to attenuate the star light. Two laser light
sources (laser 1 and laser 2) are used to simulate the starlight and the
planet light, respectively, with the optical configuration of intensity
calibration shown in Figure 4. We need to know the contrast of the two beams.
Therefore, before running the optimization code, we firstly calibrate the
intensity of the starlight and planet light sources. The step transmission
filters must be removed. A ND filter with 1$/$1000000 intensity attenuation is
inserted in the star beam. In two individual exposures, we adjust the camera
exposure time to make the light intensity in the planet and the star equal
each other. Figure 5 shows the images in two different exposure times.
According to the exposure times and ND attenuation, the contrast between the
two beams should be 85$/$0.5$/$1000000$=$1.7$\times$$10^{8}$. Once the
contrast of the two sources is known, it can be used to calculate a higher
contrast of the coronagraph. When the intensity calibration is accomplished,
we run the SPGD optimization code.

###  3.2 Optical layout of the whole system

Figure 6 shows the layout of the coronagraph optical system under test. The
starlight and planet light are simulated by two point lights using two laser
light sources ($\lambda=632.8nm$), who’s contrast are calibrated in advanced.
When we run the SPGD optimization code, the laser used to simulate the planet
light will be blocked. The starlight is suppressed when the light passes
through the pupil stop and two orthogonal step transmission filters

Due to the small number of DM actuators constraints, it is not possible to
create a large dark hole. In addition, A DM with large actuators is expensive.
Therefore, we use a LCA to replace a DM in our current test. In this test,
only one quadrant is used for high contrast imaging. A focal mask composed of
a square aperture, is located on the first focal plane (before lens L1), which
blocks three quadrants of the focal plane image and only allows one quadrant
image to go through the optical system for high-contrast imaging.

The LCA is located on a plane conjugate to the pupil. A polarizer should be
placed before the LCA in order to control incident light polarization state,
and the LCA is working as a phase-only modulator when the input light source
is linearly polarized in its vertical axis. The phase of the incoming linear
polarized light can be modulated by the LCA, as the light passes through it.
Finally, the light is imaged by a lens (L2), and a fast CCD camera takes the
PSF images.

The PSF images are taken by a 14-bit fast CCD camera and saved in computer in
real-time. Then using PSF image the SPGD algorithm evaluates the energy in the
target region of interest. According to the evaluation, the SPGD will control
the voltage of each LCA effective pixel to provide an overall phase to
minimize the energy in the target field of view. Through multiple iterations,
the expected dark-hole can be generated. Finally the computer saves the LCA
optimal voltages, which corresponds the associated phase to generate an extra
contrast improvement. And the dark hole can be locked stably in the lab
environment for at least 6 hours, which is long enough for Jupiter like
exoplanets imaging.

Figure 7: Aberrated and corrected PSFs in our testbed. The left panel shows
the PSF without focal plane mask, while the central panel shows the PSF with
focal mask. The focal plane image after the dark-hole correction is presented
in the right panel. In the right panel, the planet image (indicated by a red
circular rigion) stand out from the speckle noise in the red rectangle target
region defided by the dark hole, with a region of
$X=Y=4\lambda/D-21\lambda/D$. This monochromatic experiment is done at the
wavelength of 632.8nm.

###  3.3 Dark hole creation

The dark hole obtained using our SPGD optimization algorithm is shown in
Figure7. A dark-hole (between $4\lambda/D$ and $21\lambda/D$) is created in a
quadrant on the focal plane that is larger than JPL’s (between $4\lambda/D$
and $10\lambda/D$) and Princeton’s ($X=7-10\lambda/D$ and $Y=-3-3\lambda/D$)
(Pueyo et al. 2009). After 1000 iterations of the phase correction, the
speckle noise in the original PSF has been efficiently removed. Once the dark-
hole is created, we turn on the laser light and send the optimized LCA
phase.Then the planet image stands out the speckle noise in the dark hole
region, and can be seen clearly, as shown in the right panel of Figure 7. In
fact the intensity of the planet is about 3.5 times brighter than residual
speckle noise in the dark-hole. Therefore, the contrast in the dark hole
increases from 1$\times$$10^{-7}$ to 1.68$\times$$10^{-9}$ at $4\lambda/D$,
which demonstrates an extra gain in contrast of two orders of magnitudes
through the contribution of the dark hole using the LCA.

##  4 CONCLUSIONS AND FUTIRE WORK

In this recent test, we have presented the test results of our new 51-step
transmission filters in our high-contrast coronagraph. According to the test,
the contrast can reach $10^{-7}$ at $4\lambda/D$, without phase correction.
Then we used the LCA to generate a dark-hole based on our SPGD algorithm with
the new filters. In the dark-hole region, an extra contrast $\thicksim$100
times has been achieved. Therefore, our high-contrast coronagraph with the
step transmission filters achieves a contrast of $10^{-9}$ in the dark-hole
region. Our approach has following advantages: the wave-front control
algorithm is simple and is done based on the focal plane PSF only; the dark-
hole region created is larger, which does not need high density actuator DMs.

This paper demonstrates that our coronagraph is a promising technique for
space-based imaging earth-like planets. For future works, we will use a DM
($32\times 32$ or $64\times 64$ actuators) whose precision ($\lambda/10000$)
for phase correction is better than that of the LCA ($\lambda/100$), and the
experiment will be carried out in a space-like environment inside a vibration-
isolated and vacuum chamber. Based on the two points, the contrast should be
able to reach or better than $1\times 10^{-10}$. Further progresses will be
discussed in our future publications.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(Grant Nos. 11220101001, 11373005 and 11328302) and the ”Strategic Priority
Research Program” of the Chinese Academy of Sciences (Grant No. XDA04070600
and XDA04075200). This work was made possible through the support of a grant
from the John Templeton Foundation and National Astronomical Observatories of
Chinese Academy of Sciences. Part of the work described in this paper was
carried out at California State University Northridge, with support from the
National Science Foundation under Grant ATM-0841440.

## References

  * Borucki et al. (2013) Borucki, W. J., Agol, E., Fressin, F., et al. 2013, Science, 340, 587
  * Dong et al. (2011) Dong, B., Ren, D.-Q., & Zhang, X. 2011, RAA(Research in Astronomy and Astrophysics), 11, 997
  * Dou et al. (2009) Dou, J., Ren, D., Zhu, Y., & Zhang, X. 2009, in SPIE Optical Engineering+ Applications, 744019–744019 (International Society for Optics and Photonics)
  * Dou et al. (2010) Dou, J.-P., Ren, D.-Q., & Zhu, Y.-T. 2010, RAA(Research in Astronomy and Astrophysics), 10, 189
  * Kay et al. (2009) Kay, J. D., Pueyo, L. A., & Kasdin, N. J. 2009, in SPIE MOEMS-MEMS: Micro-and Nanofabrication, 72090G–72090G (International Society for Optics and Photonics)
  * Pepe et al. (2013) Pepe, F., Cameron, A. C., Latham, D. W., et al. 2013, Nature, 503, 377
  * Pueyo et al. (2009) Pueyo, L., Kay, J., Kasdin, N. J., et al. 2009, Appl. Opt, 48, 6296
  * Ren et al. (2012) Ren, D., Dong, B., Zhu, Y., & Christian, D. J. 2012, PASP, 124, 247
  * Ren et al. (2010) Ren, D., Dou, J., & Zhu, Y. 2010, PASP, 122, 590
  * Ren & Zhu (2007) Ren, D., & Zhu, Y. 2007, PASP, 119, 1063
  * Trauger & Traub (2007) Trauger, J. T., & Traub, W. A. 2007, Nature, 446, 771