# Optimal criteria and their asymptotic form for data selection in data-driven
reduced-order modeling with Gaussian process regression

Themistoklis P. Sapsis and Antoine Blanchard
Department of Mechanical Engineering,
Massachusetts Institute of Technology,
77 Massachusetts Ave., Cambridge, MA 02139 Corresponding author:
<EMAIL_ADDRESS>Tel: (617) 324-7508, Fax: (617) 253-8689

###### Abstract

We derive criteria for the selection of datapoints used for data-driven
reduced-order modeling and other areas of supervised learning based on
Gaussian process regression (GPR). While this is a well-studied area in the
fields of active learning and optimal experimental design, most criteria in
the literature are empirical. Here we introduce an optimality condition for
the selection of a new input defined as the minimizer of the distance between
the approximated output probability density function (pdf) of the reduced-
order model and the exact one. Given that the exact pdf is unknown, we define
the selection criterion as the supremum over the unit sphere of the native
Hilbert space for the GPR. The resulting selection criterion, however, has a
form that is difficult to compute. We combine results from GPR theory and
asymptotic analysis to derive a computable form of the defined optimality
criterion that is valid in the limit of small predictive variance. The derived
asymptotic form of the selection criterion leads to convergence of the GPR
model that guarantees a balanced distribution of data resources between
probable and large-deviation outputs, resulting in an effective way for
sampling towards data-driven reduced-order modeling.

Keywords: Optimal experimental design; Data-driven modeling; Bayesian
regression; Optimal sampling; Active learning

##  1 Introduction

Reduced-order modeling has been a cornerstone of modern computational methods.
The effectiveness of the reduced-order models relies both on their design, so
they can capture the complexity of the underlying process, and also on the
information they rely on. This fundamental information can have the form of i)
governing equations, which typically carry assumptions or simplifications of
their own, and ii) data, which is a more reliable but expensive source. The
present work involves the development of criteria for the most effective
selection of data—or associated experiments to generate this data—in order to
perform data-driven reduced-order modeling.

The literature related to data-driven reduced-order modeling is vast and spans
a great number of engineering and scientific fields ranging from fluid
mechanics [6, 12, 9, 8, 14] to structural mechanics [13, 11, 16, 1]. In the
majority of these works the assumption is plentiful data. While for many
applications this is indeed the case, there are several important problems
where plentiful data is not available, either because the associated physical
or numerical experiments are too expensive or because of the nature of the
problem, e.g., extreme events that occur rarely [20].

The scientific field that aims to tackle this issue is active learning [18, 7,
10]. A critical issue in active learning is the choice of acquisition
function, i.e., the criterion used to select which sample to query next in an
optimal manner. Acquisition functions come in various shapes and forms [7,
21], but many popular criteria suffer from severe limitations, including high
computational cost, intractability in high dimensions, and inability to
discriminate between active and idle input variables [19]. Recently, a new
class of active learning criteria was introduced, designed to take into
account the effect of the output on the selection of samples [19, 5]. This new
class of output-weighted criteria has led to significantly improved
performance in a variety of problems involving uncertainty quantification,
Bayesian optimization [3, 2], and decision making [4, 25]. However, there is
no sound theoretical understanding of its favorable properties.

In this work we rigorously show the optimal properties of output-weighted
acquisition functions when it comes to the problem of selecting samples for
data-driven reduced-order modeling in the Bayesian context. Specifically, we
employ Gaussian process regression (GPR) as the building block to develop
Bayesian surrogates or stochastic reduced-order models. Within the GPR
framework, we derive asymptotically optimal acquisition functions, where
optimality is defined in the sense of fastest convergence of the probability
density function (pdf) describing the quantity of interest. We demonstrate the
derived optimal acquisition functions in a mechanical oscillator subjected to
high-dimensional stochastic forcing resulting in non-Gaussian heavy tails, as
well as the reduced-order modeling of a beam under axial and transverse
stochastic loads that result in buckling and bending.

##  2 Problem setup

Our aim is to build a data-driven reduced-order model that will capture the
behavior of an output quantity $y\in\mathbb{R}$ (assumed scalar for
simplicity) with respect to an input variable $x\in\mathbb{R}^{n}$. The input
can represent, for example, initial conditions or parameters governing the
evolution of a dynamical system, while the output is any quantity of interest
that depends on the input variable. We also assume that the input variable has
a prescribed probability distribution function, $p_{x}({x})$.

We are given a set of input datapoints $X={\left\\{x_{i}\right\\}}_{i=1}^{N}$
and corresponding output values

$Y=[y(x_{1}),...,y(x_{N})]^{\mathsf{T}},$

and the objective is to identify a criterion that allow us to select the next
input point $x_{N+1}\triangleq h$ that will supplement—together with the
corresponding output $y(x_{N+1})$—the dataset $\mathcal{D}=\\{X,Y\\}$. There
is a plethora of criteria in the active-learning field (see, e.g., [7] for a
review), the majority of which are designed to select the next input $h$ that
either minimizes the uncertainty of the reduced-order model, i.e., the
posterior variance $\sigma^{2}_{y}(x|\mathcal{D})$, or maximizes the
information content between input and output variables. However, neither of
these methods gives appropriate attention on increasing the accuracy of the
resulting model in the regions of the input space where it is needed the most:
inputs that result in large deviations of the output variable $y$ from its
expected value.

To accomplish this goal we will focus on accelerating the convergence of the
model output pdf, $p_{y}(s|\mathcal{D})$. In the context of reduced-order
modeling this is the appropriate quantity to consider as it encodes the
probabilistic information about the large deviations of $y$, while the input
variables have been marginalized. Therefore, while criteria based on the
posterior variance aim to minimize the error of the model for all or certain
input variables, a criterion based on the output pdf explicitly targets
regions of the input space that contribute to the output pdf. Moreover, to
emphasize regions associated with large deviations of the output variable $y$,
we will define convergence in terms of the logarithm of
$p_{y}(s|\mathcal{D})$, using the below function to measure the discrepancy
between two pdfs $p_{1}$ and $p_{2}$:

$\mathbb{D}(p_{1},p_{2})=\int_{S_{y}}\left|\log p_{1}({s})-\log
p_{2}({s})\right|\mathrm{d}{s}.$ (1)

In the above, $S_{y}$ is a finite domain of $y$ (the output variable), over
which we aim to build the reduced-order model. Note that this is different
from the Kullback–Leibler divergence, which does not give the same emphasis to
low probability events (associated with large deviations of $y$). Moreover, it
can be easily shown that the above function defines a metric. In what follows,
we first provide a quick review of GPR and its basic properties, and
subsequently formulate the optimality condition that our sampling criteria
will satisfy.

###  2.1 Review of Gaussian process regression

We use GPR to build a surrogate for the unknown function
$y(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$. The idea is to utilize a Gaussian
prior on the function $y$, i.e.,

$y_{0}\sim\text{GP}(m(x),k(x,x^{\prime})),$ (2)

with known mean $m(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$ and covariance
function $k(x,x^{\prime}):\mathbb{R}^{n\times n}\rightarrow\mathbb{R}$ assumed
to be positive-definite. We typically set the mean function $m$ to be
identically zero. As for the covariance kernel $k$, we often use the Gaussian
covariance:

$k(x,x^{\prime})=\sigma_{k}^{2}\exp\left(-\frac{\left\|x-x^{\prime}\right\|^{2}}{2\lambda^{2}}\right).$
(3)

Conditioning the Gaussian process on the available datapoints, we obtain the
predictive process

$\displaystyle y_{*}\sim\text{GP}(\bar{y}(x),\bar{k}(x,x^{\prime})),$ (4)

where the predictive mean, $\bar{y}(x)$ and predictive covariance
$\bar{k}(x,x^{\prime})$ are given in terms of the datapoints:

$\displaystyle\begin{split}\bar{y}(x)&=k(x,X)K(X,X)^{-1}Y\\\
\bar{k}(x,x^{\prime})&=k(x,x^{\prime})-k(x,X)K(X,X)^{-1}k(X,x^{\prime}),\end{split}$
(5)

where $k(x,X)=[k(x,x_{1}),...,k(x,x_{N})]\in\mathbb{R}^{N}$ and
$K(X,X)=\left\\{k(x_{i},x_{j})\right\\}_{i,j=1}^{N}\in\mathbb{R}^{N\times N}$.

The native space for GPR schemes is the reproducing kernel Hilbert space
corresponding to the kernel $k$, defined as follows:

######  Definition 2.1.

A Hilbert space $\mathcal{H}_{k}$ of functions
$y:\mathbb{R}^{n}\rightarrow\mathbb{R}$, with inner product
$\left\langle\cdot,\cdot\right\rangle_{\mathcal{H}_{k}}$, is called the
reproducing kernel Hilbert space (RKHS) corresponding to a symmetric,
positive-definite kernel $k$ if

  1. 1.

for all $x\in\mathbb{R}^{n}$, $k(x,x^{\prime})$, as a function of its second
argument, $x^{\prime}$, belongs to $\mathcal{H}_{k}$; and

  2. 2.

for all $x\in\mathbb{R}^{n}$ and $y\in\mathcal{H}_{k}$, $\left\langle
y,k(x,\cdot)\right\rangle_{\mathcal{H}_{k}}=y(x)$.

For special choices of kernels $k$, the RKHS can be characterized through its
spectrum. Specifically, we have the following theorem for the characterization
of the RKHS in the case of Gaussian kernels:

######  Theorem 2.2 (Wendland, 2004, Theorem 10.12, [24]).

Let
$k(x,x^{\prime})=\sigma_{k}^{2}\exp(-{\left\|x-x^{\prime}\right\|^{2}}/{2\lambda^{2}})$
be the squared-exponential kernel. The corresponding RKHS $\mathcal{H}_{k}$
can be written as

$\displaystyle\mathcal{H}_{k}=\left\\{y\in L_{2}(\mathbb{R}^{n})\cap
C(\mathbb{R}^{n}):\left\|y\right\|_{\mathcal{H}_{k}}=\frac{1}{c_{0}}\int|\mathcal{F}[y](\omega)|^{2}\exp(\lambda^{2}\left\|\omega\right\|^{2}/2)\,\mathrm{d}\omega<\infty\right\\},$
(6)

where $c_{0}$ is a constant that depends on $n$ and $\lambda$, and
$\mathcal{F}$ is the Fourier transform.

The above results shows that for any $y\in\mathcal{H}_{k}$, the magnitude of
its Fourier transform $|\mathcal{F}[y](\omega)|$ decays exponential fast as
$|\omega|\rightarrow\infty$ and the speed of decay increases with $\lambda$.
Analogous results exist for the case of Matérn kernels [24].

A fundamental property of GPR schemes is the fact that one can obtain a priori
estimates for the accuracy of the surrogate model. In particular, for the case
of GPR schemes, we have the following error estimate:

######  Proposition 2.3 (Stuart & Teckentrup, 2018, Proposition 3.5, [22]).

Suppose that $\bar{y}(x)$ and $\bar{k}(x,x^{\prime})$ are given by the GPR
scheme (5). Then

$\displaystyle\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}|y(x)-\bar{y}(x)|=\bar{k}(x,x)^{\frac{1}{2}}\triangleq\bar{\sigma}(x),$
(7)

where the supremum occurs when the functions $y(\cdot)$ and $\bar{k}(\cdot,x)$
are linearly dependent.

Based on this result which involves the accuracy of the surrogate map, we will
derive the corresponding results for the pdf of the output variable,
$p_{y}(s|\mathcal{D})$.

###  2.2 Optimality condition for data selection

We formulate an active-sampling criterion that aims directly for the
convergence of the output pdf, $p_{y}(s|\mathcal{D})$. Let
$h\in\mathbb{R}^{n}$ be the new candidate input point. As output, we employ
the approximation by the surrogate model, $\bar{y}(h)$. In this way, we have
the augmented dataset $\mathcal{D^{\prime}}=\\{[X,h],[Y,\bar{y}(h)]\\}$, which
results in the same predictive mean, $\bar{y}(x)$, and a new predictive
covariance $\bar{k}^{\prime}(x,x^{\prime};h)$. Ideally, we would want the new
input to be chosen by minimizing the distance

$\mathbb{D}(p_{y},p_{\bar{y}^{\prime}})=\int_{S_{y}}\left|\log p_{y}({s})-\log
p_{\bar{y}^{\prime}}({s}|\mathcal{D},h)\right|\mathrm{d}{s}.$ (8)

However, this is not possible as $p_{y}$ is a priori unknown. To this end, we
will use as a selection criterion the supremum of the above distance over the
unit sphere of the functional space $y\in\mathcal{H}_{k}$ (we fix the norm of
the unknown function without loss of generality). Therefore the selection
criterion takes the form of minimizing the acquisition function

$Q(h|\mathcal{D})\triangleq\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}\mathbb{D}(p_{y},p_{\bar{y}^{\prime}})=\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}\int_{S_{y}}\left|\log
p_{y}({s})-\log p_{\bar{y}^{\prime}}({s}|\mathcal{D},h)\right|\mathrm{d}{s}.$
(9)

While the criterion is targeting the output pdf, it is not easily computable,
especially in high dimensions. The rest of the paper aims to derive a
computable version appropriate for high-dimensional input spaces.

##  3 Asymptotically optimal criterion for data selection

Our efforts focus on obtaining a computable version of the criterion (9). We
plan to achieve this by assuming small variance
$\bar{\sigma}^{2}(x)=\bar{k}(x,x)$. We first recall an asymptotic result that
connects the error for a map and the error between the induced pdfs defined by
the corresponding maps.

######  Theorem 3.1 (Mohamad & Sapsis, 2018, Theorem 2, [15]).

Let $\hat{y}(x)$ and $y(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$ be two
continuous functions with difference $\Delta y({x})=\hat{y}(x)-y(x)$, which is
assumed to be small. Let also $p_{x}({x})$ be the probability density function
of the random vector $x\in\mathbb{R}^{n}$. The difference between the induced
pdfs for $\hat{y}$ and $y$ has the following asymptotic behavior:

$p_{y}({s})-p_{\hat{y}}({s})=-\frac{\mathrm{d}}{\mathrm{d}{s}}\int\displaylimits_{{\hat{y}}({x})={s}}p_{x}({x})\Delta
y({x})\,\mathrm{d}{x}+\mathcal{O}(|\Delta y|^{2}).$

Building on this result, we have the main theorem that characterizes the
supremum in the RKHS between the pdf induced by the surrogate approximation
and the pdf induced by the exact map:

######  Theorem 3.2.

Let $y(x)\in\mathcal{H}_{k}$ be an arbitrary function and its GPR
approximation with kernel $k$ given by ${\bar{y}(x)}$ with corresponding
variance $\bar{\sigma}^{2}(x)$. Then the following property holds for small
$\bar{\sigma}$:

$\displaystyle\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}\int_{S_{y}}\left|\log
p_{\bar{y}}({s})-\log
p_{y}({s})\right|\mathrm{d}{s}=\int\displaylimits_{\bar{y}^{-1}(S_{y})}\frac{p_{x}({x})|p^{\prime}_{y}({\bar{y}}({x}))|}{p_{y}({\bar{y}}({x}))^{2}}\bar{\sigma}({x})\,\mathrm{d}{x}+\mathcal{O}(\bar{\sigma}^{2}).$
(10)

Proof:  Utilizing Theorem 3.1, we have for the case where ${\bar{y}}$ is close
to $y$ (i.e., $\Delta y/y=(\hat{y}-y)/y\ll 1$)

$\displaystyle\log p_{\bar{y}}({s})-\log p_{y}({s})$
$\displaystyle=\frac{p_{\bar{y}}({s})-p_{y}({s})}{p_{y}({s})}+\mathcal{O}(|\Delta
p_{y}|^{2})$
$\displaystyle=-\frac{\frac{\mathrm{d}}{\mathrm{d}{s}}\int\displaylimits_{{\bar{y}}({x})={s}}p_{x}({x})\Delta
y({x})\,\mathrm{d}{x}}{p_{y}({s})}+\mathcal{O}(|\Delta y|^{2}).$

Expressing the right-hand side as a volume integral with a delta function and
using properties of generalized derivatives, we have

$\displaystyle\frac{\frac{\mathrm{d}}{\mathrm{d}{s}}\int\displaylimits_{{\bar{y}}({x})={s}}p_{x}({x})\Delta
y({x})\,\mathrm{d}{x}}{p_{y}({s})}$
$\displaystyle=\frac{\frac{\mathrm{d}}{\mathrm{d}{s}}\int\displaylimits
p_{x}({x})\Delta y({x})\delta({s-{\bar{y}}({x})})\,\mathrm{d}{x}}{p_{y}({s})}$
$\displaystyle=\int\displaylimits\frac{p_{x}({x})\Delta
y({x})\delta^{\prime}({s-{\bar{y}}({x})})}{p_{y}({s})}\,\mathrm{d}{x}$
$\displaystyle=\int\displaylimits\frac{p_{x}({x})\Delta
y({x})p^{\prime}_{y}({s})\delta({s-{\bar{y}}({x})})}{p^{2}_{y}({s})}\,\mathrm{d}{x}$
$\displaystyle=p^{\prime}_{y}({s})\int\displaylimits\frac{p_{x}({x})\Delta
y({x})\delta({s-{\bar{y}}({x})})}{p^{2}_{y}({s})}\,\mathrm{d}{x}.$

We note that $p_{x}({x})\delta({s-{\bar{y}}({x})})/p^{2}_{y}({s})\geq 0$. By
employing Proposition 2.3 we obtain the tight bound

$\displaystyle\int\frac{p_{x}({x})\Delta
y({x})\delta({s-{\bar{y}}({x})})}{p^{2}_{y}({s})}\,\mathrm{d}{x}\leq\int\displaylimits\frac{p_{x}({x})\bar{\sigma}({x})\delta({s-{\bar{y}}({x})})}{p^{2}_{y}({s})}\,\mathrm{d}{x},$

where equality holds for the case where the functions $y(\cdot)$ and
$\bar{k}(\cdot,x)$ are linearly dependent [22].

Combining the above, we have

$\displaystyle\left|\log p_{\bar{y}}({s})-\log p_{y}({s})\right|$
$\displaystyle\leq|p^{\prime}_{y}({s})|\int\displaylimits\frac{p_{x}({x})\bar{\sigma}({x})\delta({s-{\bar{y}}({x})})}{p^{2}_{y}({s})}\,\mathrm{d}{x}+\mathcal{O}(\bar{\sigma}^{2}).$

Integrating over ${s}$, we have the final result:

$\displaystyle\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}\int_{S_{y}}\left|\log
p_{\bar{y}}({s})-\log p_{y}({s})\right|\mathrm{d}s$
$\displaystyle=\int_{S_{y}}|p^{\prime}_{y}({s})|\int\displaylimits\frac{p_{x}({x})\bar{\sigma}({x})\delta({s-{\bar{y}}({x})})}{p^{2}_{y}({s})}\,\mathrm{d}{x}\,\mathrm{d}{s}+\mathcal{O}(\bar{\sigma}^{2})$
$\displaystyle=\int\displaylimits_{\bar{y}^{-1}(S_{y})}\frac{p_{x}({x})|p^{\prime}_{y}({{\bar{y}}({x})})|}{p^{2}_{y}({\bar{y}}({x}))}\bar{\sigma}({x})\,\mathrm{d}{x}+\mathcal{O}(\bar{\sigma}^{2}).$

This completes the proof. $\blacksquare$

Next, we utilize this asymptotic form to reformulate the data selection
criterion (9). Specifically, we apply the above theorem for the augmented
dataset $\mathcal{D}^{\prime}$ and obtain the asymptotic reformulation of the
selection criterion:

$\displaystyle\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}\int_{S_{y}}\left|\log
p_{y}({s})-\log p_{\bar{y}^{\prime}}({s}|\mathcal{D},h)\right|\mathrm{d}{s}$
$\displaystyle=\int\displaylimits_{\bar{y}^{-1}(S_{y})}\frac{p_{x}({x})|p^{\prime}_{y}({\bar{y}}({x}))|}{p_{y}({\bar{y}}({x}))^{2}}\bar{\sigma}({x;h})\,\mathrm{d}{x}+\mathcal{O}(\sigma^{2}),$
(11)

where $\bar{\sigma}({x;h})=\bar{k}(x,x;h)$ is the predictive variance based on
the augmented dataset $\mathcal{D}^{\prime}$. We note that the right-hand side
involves the pdf $p_{y}$, which is unknown but can always be approximated by
$p_{\bar{y}}$ with negligible error, given the assumptions of Theorem 3.2.
This gives us the final asymptotic approximation:

$\displaystyle
Q(h|\mathcal{D})=\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}\int_{S_{y}}\left|\log
p_{y}({s})-\log p_{\bar{y}^{\prime}}({s}|\mathcal{D},h)\right|\mathrm{d}{s}$
$\displaystyle\simeq\int\displaylimits_{\bar{y}^{-1}(S_{y})}\frac{p_{x}({x})|p^{\prime}_{\bar{y}}({\bar{y}}({x}))|}{p_{\bar{y}}({\bar{y}}({x}))^{2}}\bar{\sigma}({x;h})\,\mathrm{d}x.$
(12)

It is worth emphasizing the term in the denominator, which promotes sampling
of regions associated with low probability, i.e., large deviations of the
output of the reduced-order model.

We can simplify the right-hand side further by using the Cauchy–Schwarz
inequality (with weight $\frac{p_{x}({x})}{p_{\bar{y}}({\bar{y}}({x}))}$) to
obtain the less conservative upper bound:

$Q(h|\mathcal{D})\simeq\int\displaylimits_{\bar{y}^{-1}(S_{y})}\frac{p_{x}({x})|p^{\prime}_{\bar{y}}({{\bar{y}}({x})})|}{p^{2}_{\bar{y}}({\bar{y}}({x}))}\bar{\sigma}({x;h})\,\mathrm{d}{x}\leq
c\left[\int\displaylimits_{\bar{y}^{-1}(S_{y})}\frac{p_{x}({x})}{p_{\bar{y}}({\bar{y}}({x}))}\bar{\sigma}^{2}({x};h)\,\mathrm{d}x\right]^{\frac{1}{2}},$
(13)

where

$c=\left[\int\displaylimits_{\bar{y}^{-1}(S_{y})}\frac{p_{x}({x})p^{\prime
2}_{\bar{y}}({{\bar{y}}({x})})}{p^{3}_{\bar{y}}({\bar{y}}({x}))}\,\mathrm{d}{x}\right]^{\frac{1}{2}}$
(14)

is a constant that depends only $p_{x}(x)$, $p_{\bar{y}}(y)$, and
$\bar{y}(x)$, i.e., not on $h$. This form, also referred as output-weighted
(or likelihood-weighted) criterion, is appropriate for computations involving
even high-dimensional input spaces, since it allows for the analytical
computation of $\sigma^{2}({x};h)$ in terms of simpler integrals [5]. It has
been studied numerically in recent papers [19, 5, 3], showing significantly
favorable convergence properties compared with existing active-learning
criteria.

###  3.1 The case of extreme-event quantiles

For a wide range of applications, the focus is on characterizing the
probability of exceeding a certain level, $P_{y}(s_{*})=P[y\leq s_{*}]$,
rather than characterizing the full pdf $p_{y}(s)$ (the case of non-exceeding
probability can be addressed in a similar fashion). For this type of problems
the appropriate selection criterion is minimizing the acquisition function

$R(h|\mathcal{D},s_{*})\triangleq\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}|P_{\bar{y}^{\prime}}({s_{*}}|\mathcal{D},h)-P_{y}({s_{*}})|.$
(15)

For this case, we have the following asymptotic result for small $\sigma$:

######  Theorem 3.3.

Let $y(x)\in\mathcal{H}_{k}$ and its GPR approximation with kernel $k$ given
by ${\bar{y}(x)}$ with variance $\bar{\sigma}^{2}(x)$. Then, for any given
$s_{*}$, the following property holds for small $\bar{\sigma}$:

$\displaystyle\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}|P_{\bar{y}}({s_{*}})-P_{y}({s_{*}})|=\int\displaylimits_{\bar{y}({x})={s_{*}}}\bar{\sigma}({x})p_{x}({x})\,\mathrm{d}x+\mathcal{O}(\bar{\sigma}^{2}).$

Proof:  The starting point is Theorem 3.1. We integrate from $-\infty$ to
$s_{*}$ to obtain the cumulative distribution function on the left-hand side,
when ${\bar{y}}$ is close to $y$ (i.e., $\Delta y/y=(\hat{y}-y)/y\ll 1$):

$\displaystyle
P_{\bar{y}}({s})-P_{y}({s})=-\int\displaylimits_{{\bar{y}}({x})={s}}p_{x}({x})\Delta
y({x})\,\mathrm{d}{x}+\mathcal{O}(|\Delta y|^{2}).$

We then employ Proposition 2.3 and bound the difference $\Delta y$ on the
right-hand side. This completes the proof. $\blacksquare$

Applying the above theorem to the augmented dataset $\mathcal{D}^{\prime}$, we
obtain the asymptotic form for the optimal selection criterion involving
quantiles:

$\displaystyle
R(h|\mathcal{D},s_{*})\triangleq\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}|P_{\bar{y}^{\prime}}({s_{*}}|\mathcal{D},h)-P_{y}({s_{*}})|=\int\displaylimits_{\bar{y}({x})={s_{*}}}\bar{\sigma}({x;h})p_{x}({x})\,\mathrm{d}x+\mathcal{O}(\bar{\sigma}^{2}).$
(16)

We observe that for this case, the optimal data selection criterion takes a
different form, focusing primarily on reducing the error around the contour
$\bar{y}(x)=s_{*}$. This is not surprising given that the extreme-event
quantile is essentially equivalent to a classification problem, so what is
most important is to have low error around the a priori unknown contour of
interest.

###  3.2 Convergence of spatial approximation error

Given an active-learning criterion, it is important to characterize the
convergence of the spatial error of the GPR approximation. This can be
obtained with the help of Proposition 2.3 and standard results of measurable
functions. Here we give a proof for the convergence properties using the
sampling criterion appearing on the right-hand side in (13). Specifically, we
have the following result:

######  Theorem 3.4.

Suppose that ${\bar{y}_{N}}$ and $\bar{\sigma}_{N}$ are given by a GPR with
kernel $k$, using $N$ samples. Moreover, assume that sampling is performed
using the optimal criterion so that

$\lim_{N\rightarrow\infty}\int\frac{p_{x}({x})}{p_{\bar{y}_{N}}({\bar{y}_{N}}({x}))}{\bar{\sigma}_{{N}}^{2}}({x})\,\mathrm{d}{x}=0.$

Then we have convergence in measure, i.e.,

$\lim_{N\rightarrow\infty}P\left[x:\sup_{\left\|y\right\|_{\mathcal{H}_{k}}=1}({\bar{y}_{N}}({x})-{y}({x}))^{2}\leq
q\frac{p_{\bar{y}_{N}}({\bar{y}_{N}}(x))}{p_{x}({x})}\right]=1,$

for every $q>0$.

Proof: From standard results of measurable functions [23], for every sequence
of functions $\phi_{k}(x)\in L^{p}$ ($1\leq p\leq\infty$) for which
$\|\phi_{k}(x)\|_{p}\rightarrow 0$, we have convergence in measure:

$\displaystyle\lim_{k\rightarrow\infty}P[x:|\phi_{k}(x)|\geq q]=0,$ (17)

for every $q>0$. Applying this result to the sequence
$p_{x}({x})\bar{\sigma}_{{N}}^{2}({x})/p_{\bar{y}_{N}}({\bar{y}_{N}}({x}))$,
we obtain

$\displaystyle\lim_{N\rightarrow\infty}P\left[x:\frac{p_{x}({x})}{p_{\bar{y}_{N}}({\bar{y}_{N}}({x}))}{\bar{\sigma}_{{N}}^{2}}({x})\leq
q\right]=1.$

We then utilize Proposition 2.3, which immediately leads to the desired
result. $\blacksquare$

This result provides a description of the spatial convergence properties for
the approximation error of the reduced-order model. Specifically, it shows
that the convergence is accelerated in regions of the input space that are i)
most probable according to the pdf $p_{x}(x)$, and ii) associated with small
probability of the output pdf $p_{\bar{y}_{N}}({\bar{y}_{N}}(x))$, i.e., large
deviations of the output $y$. This type of convergence guarantees a balanced
distribution of resources between the most probable inputs and those that
result in large deviations for the output, resulting in an effective way of
sampling towards data-driven reduced-order modeling.

##  4 Numerical illustration

We demonstrate the optimal sampling criteria in a mechanical oscillator
subject to stochastic forcing and in the reduced-order modeling of a beam
under axial and transverse stochastic loads. We consider the optimal sampling
criteria appearing on each side of the inequality in (13). The left-hand side
will be referred to as the “B” criterion, and the right-hand side as “IVR-LW”
as it is strictly equivalent to the eponym criterion introduced in [19, 5]. We
also consider two criteria commonly used in the literature which do not
account for the importance of the output relative to the input; namely,
uncertainty sampling (US) and input-weighted integrated variance reduction
(IVR-IW), whose definitions can be found in [19, 5]. For each example below,
we run 100 Bayesian experiments, each differing in the choice of the $n+1$
points making up the initial dataset. Observations are assumed to be corrupted
by Gaussian noise with zero mean and unknown (i.e., to be learned) variance
$\sigma_{n}^{2}$. Performance at each iteration is evaluated using the median
of $\mathbb{D}(p_{y},p_{\bar{y}_{N}})$ across the 100 randomized experiments.
Note that we employ this particular quantity to measure performance so we can
better emphasize the accuracy in the tails of the resulting probability
distributions.

###  4.1 Forced nonlinear oscillator exhibiting extreme events

We begin with the stochastic oscillator of Mohamad and Sapsis [15],

$\ddot{u}+\delta\dot{u}+F(u)=\xi(t),\quad t\in[0,T],$ (18)

where $u(t)\in\mathbb{R}$ is the state variable, $F$ a nonlinear restoring
force, and $\xi(t)$ a stationary stochastic process which we parametrize using
a Karhunen–Loève expansion with $n$ modes:

$\xi(t)\simeq x\Phi(t),\ $

where $\\{\Lambda,\Phi(t)\\}$ contains the first $n$ eigenpairs of the
correlation matrix, and $x\in\mathbb{R}^{n}$ is a vector of random
coefficients having mean zero and diagonal covariance matrix $\Lambda$. The
quantity of interest is taken to be the mean value of $u(t)$ over the interval
$[0,T]$. (Further details about system parameters can be found in [15, 5]).

We consider the case $n=2$ as it allows visual comparison of the decisions
made by the sampling criteria as more points are being acquired. Even with
$n=2$, the output pdf has heavy tails (see figure 1), a consequence of the
strong nonlinearity in (18). For $\sigma_{n}^{2}=10^{-3}$, figure 1 shows that
the derived optimal criteria accelerate convergence of the output pdf quite
dramatically. The error for B remains close to, but always slightly below,
that for IVR-LW, which is consistent with the mathematical derivation laid out
in the previous section.

Figure 1: For the stochastic oscillator (18) with $n=2$, contour plot of the
output $y$ (left) and its pdf (center), and performance of several sampling
criteria for $\sigma_{n}^{2}=10^{-3}$ (right). The error bands indicate one
half of the median absolute deviation.

To explain the success of the proposed optimal criteria, we investigate the
decisions made by US, B, and IVR-LW in the case where observations are
noiseless and $\sigma_{n}^{2}$ is set to zero in the GPR model (i.e., it is
not learned from data). Consistent with [19, 5], figure 2 shows that US
attempts to reduce uncertainty somewhat evenly across the space as it has no
mechanism to discriminate between relevant and irrelevant regions. With IVR-
IW, the algorithm does not explore beyond the center region where $p_{x}$ is
large, and consequently the interesting regions are not visited. On the other
hand, both B and IVR-LW decide to focus on a diagonal band, with the former
being even more surgical and localized than the latter. Specifically, IVR-LW
focuses on input regions with important probability as well as those input
regions associated with large outputs. In this way the resulting surrogates
predict the output statistics much better than with US.

(a)  (b)

(c)  (d)

Figure 2: For the stochastic oscillator (18) with $n=2$ and
$\sigma_{n}^{2}=0$, progression of the sampling algorithm for several criteria
after 60 iterations; the contours denote the posterior mean of the GPR model,
the open squares the initial dataset, and the filled circles the optimized
samples.

###  4.2 Buckling of a beam under stochastic axial and transverse excitation

Next, we consider the case of a beam of length $l$ subject to both axial and
transverse stochastic loads. The linearized equation of motion takes the form

$\displaystyle\frac{\partial^{2}w}{\partial
t^{2}}+2\zeta\omega_{0}\frac{\partial w}{\partial
t}+\omega_{0}^{2}\frac{\partial^{4}w}{\partial
x^{4}}+P(t)\frac{\partial^{2}w}{\partial x^{2}}=R(x,t),$ (19)

where $P(t)$ is random axial load and $R(x,t)$ is a random distributed
transverse load [17]. We assume that both functions have zero-mean Gaussian
statistics and prescribed spectra. We also assume that the beam has pin
boundary conditions at both ends. In this case and assuming weak damping, the
solution can be expressed as

$\displaystyle w(x,t)=\sum_{j=1}^{\infty}f_{j}(t)\sin(j\pi x/l),$ (20)

which results in a set of modal equations:

$\displaystyle\ddot{f}_{j}+2\zeta\omega_{0}\dot{f}+\omega_{j}^{2}[1-P(t)/c_{j}]f=R_{j}(t),\qquad
j=1,2,...$ (21)

where

$\displaystyle\omega_{j}^{2}=\omega_{0}^{2}(j\pi/l)^{4},\quad
c_{j}=(j\pi/l)^{4},\quad\text{and}\quad
R_{j}(t)=\frac{2}{l}\int_{0}^{l}\sin(j\pi x/l)R(x,t)\,\mathrm{d}x.$ (22)

We define as quantity of interest the maximum absolute displacement at $x=l/4$
over a prescribed time interval $[0,T]$:

$\displaystyle
y=\max_{t\in[0,T]}\left|w(l/4,t)\right|=\max_{t\in[0,T]}\left|\sum_{j=1}^{J}\sin(j\pi/4)f_{j}(t)\right|.$
(23)

We use $\sigma_{\xi}^{2}\exp[-t^{2}/(2\ell_{\xi}^{2})]$ for the correlation
function of $P(t)$ and each $R_{j}(t)$, which are expanded using
Karhunen–Loève expansions with $n_{\textit{KL}}$ modes. The search space,
therefore, has dimension $n=J(n_{\textit{KL}}+1)$. We use parameters
$\sigma_{\xi}=20$, $\ell_{\xi}=0.1$, and $T=5$.

For simplicity, we consider a modal truncation of $J=1$ and
$n_{\textit{KL}}=1$, leading to a two-dimensional search space. For the
parameters considered, figure 3 shows that the pdf of the output has a heavy
right tail, to which the B and IVR-LW criteria converge more quickly than US
and IVR-IW. The reason is that the extreme displacement values are found in
low-probability areas of the search space (i.e., small $p_{x}$), which US and
IVR-IW have no mechanism to discover.

Figure 3: For the beam under random load with $n=2$ ($J=1$,
$n_{\textit{KL}}=1$), contour plot of the output $y$ (left) and its pdf
(center), and performance of several sampling criteria for
$\sigma_{n}^{2}=10^{-3}$ (right). The error bands indicate one half of the
median absolute deviation.

##  5 Conclusions

We have derived optimal acquisition functions for active-learning schemes
utilized for reduced-order modeling based on Gaussian process regression. The
key feature of these optimal acquisition functions is a mechanism that targets
the output pdf of the surrogate model. The derivation begins by selecting each
sample so that the distance between the exact output pdf and the approximated
output pdf obtained from the reduced-order model is minimized. Given that the
exact pdf is a priori unknown, a supremum of this distance over the native
Hilbert space for the Gaussian process regression is considered. The resulting
selection criterion, although optimal, it is generally computationally
intractable. We show that this difficulty can be overcome by deriving
successive asymptotic upper bounds, resulting in a sampling criterion which
can be evaluated analytically along with its gradients. In addition to the
data-selection criterion, we derive the corresponding bound for the
approximation error of the resulted reduced-order model. This result shows
that our approach enables the reduced-order model to find the optimal trade-
off between most likely and most interesting (i.e. large-deviation) output for
the process at hand, thereby greatly enhancing the effectiveness of the
algorithm. Numerical results confirm the derived analytical findings.

Note that in this work and for the sake of simplicity we have chosen to
represent the reduced-order model in the form of a function from the parameter
space to the output space. The present framework can be adapted to represent
the reduced-order model in the form of a low-dimensional dynamical system,
focusing e.g. on capturing a specific mechanism of the full dynamics. We leave
this topic as a possible direction for future work.

#### Acknowledgments

The authors acknowledge support from the Air Force Office of Scientific
Research (MURI Grant No. FA9550-21-1-0058), the Defense Advanced Research
Projects Agency (Grant No. HR00112110002), and the Office of Naval Research
(Grant No. N00014-21-1-2357).

## References

  * [1] T. Bertalan, F. Dietrich, I. Mezić, and I. G. Kevrekidis.  On learning Hamiltonian systems from data.  Chaos, 29(12):1–9, 2019.
  * [2] A. Blanchard, G. C. Maceda, D. Fan, Y. Li, Y. Zhou, B. Noack, and T. P. Sapsis.  Bayesian optimization for active flow control.  Acta Mechanica Sinica, 2021.
  * [3] A. Blanchard and T. Sapsis.  Bayesian Optimization with Output-Weighted Importance Sampling.  Journal of Computational Physics, 425:109901, 2021.
  * [4] A. Blanchard and T. Sapsis.  Informative Path Planning for Extreme Anomaly Detection in Environment Exploration and Monitoring.  arXiv, 2021.
  * [5] A. Blanchard and T. Sapsis.  Output-Weighted Optimal Sampling for Bayesian Experimental Design and Uncertainty Quantification.  SIAM ASA J. of Uncertainty Quantification, 9:564, 2021.
  * [6] S. Brunton, B. Noack, and P. Koumoutsakos.  Machine Learning for Fluid Mechanics.  Ann. Rev. Fluid Mech., 52:477–508, 2020.
  * [7] K. Chaloner and I. Verdinelli.  Bayesian experimental design: A review.  Statistical Science, 10(3):273–304, 1995.
  * [8] D. Fernex, R. Semaan, M. Albers, P. S. Meysonnat, W. Schröder, and B. R. Noack.  Actuation response model from sparse data for wall turbulence drag reduction.  Physical Review Fluids, 2020.
  * [9] O. Ghattas and K. Willcox.  Learning physics-based models from data: perspectives from inverse problems and model reduction.  Acta Numerica, 30:445–554, 2021.
  * [10] R. B. Gramacy and H. K. Lee.  Adaptive design and analysis of supercomputer experiments.  Technometrics, 2009.
  * [11] S. Jain and G. Haller.  How to Compute Invariant Manifolds and their Reduced Dynamics in High-Dimensional Finite-Element Models?  pages 1–40, 2021.
  * [12] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3(6):422–440, 2021.
  * [13] G. Kerschen, K. Worden, A. F. Vakakis, and J. C. Golinval.  Past, present and future of nonlinear system identification in structural dynamics.  Mechanical Systems and Signal Processing, 20(3):505–592, 2006.
  * [14] Z. Ma and W. Pan.  Data-driven nonintrusive reduced order modeling for dynamical systems with moving boundaries using Gaussian process regression.  Computer Methods in Applied Mechanics and Engineering, 373:113495, 2021.
  * [15] M. A. Mohamad and T. P. Sapsis.  Sequential sampling strategy for extreme event statistics in nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 115(44):11138–11143, 2018.
  * [16] K. J. Moore, A. Mojahed, M. Kurt, M. Eriten, D. M. McFarland, L. A. Bergman, and A. F. Vakakis.  Advanced nonlinear system identification for modal interactions in nonlinear structures: A review, volume 94.  Springer International Publishing, 2019.
  * [17] A. H. Nayfeh and D. T. Mook.  Nonlinear Oscillations.  Wiley-Interscience, New York, 1984.
  * [18] J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn.  Design and analysis of computer experiments.  Statistical Science, 1989.
  * [19] T. Sapsis.  Output-weighted optimal sampling for Bayesian regression and rare event statistics using few samples.  Proceedings of the Royal Society A, 476:20190834, 2020.
  * [20] T. P. Sapsis.  Statistics of Extreme Events in Fluid Flows and Waves.  Annual Review of Fluid Mechanics, 53:85–111, 2021.
  * [21] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. De Freitas.  Taking the human out of the loop: A review of Bayesian optimization, 2016.
  * [22] A. M. Stuart and A. L. Teckentrup.  Posterior Consistency for Gaussian Process Approximations of Bayesian Posterior Distributions.  Math. Comp., 87:721–753, 2018.
  * [23] B. Z. Vulikh.  A brief course in the theory of functions of a real variable.  Mir Publishers, 1976.
  * [24] H. Wendland.  Scattered Data Approximation.  Cambridge Univ Press, 2004.
  * [25] Y. Yang, A. Blanchard, T. Sapsis, and P. Perdikaris.  Output-Weighted Sampling for Multi-Armed Bandits with Extreme Payoffs.  arXiv, 2021.


# Scattering correcting wavefront shaping for three-photon microscopy

Bernhard Rauer<EMAIL_ADDRESS>Laboratoire Kastler Brossel, ENS-
Université PSL, CNRS, Sorbonne Université, Collège de France, 24 rue Lhomond,
75005 Paris, France  Hilton B. de Aguiar  Laboratoire Kastler Brossel, ENS-
Université PSL, CNRS, Sorbonne Université, Collège de France, 24 rue Lhomond,
75005 Paris, France  Laurent Bourdieu  Institut de Biologie de l’ENS (IBENS),
École Normale Supérieure, CNRS, INSERM, Université PSL, Paris, France  Sylvain
Gigan  Laboratoire Kastler Brossel, ENS-Université PSL, CNRS, Sorbonne
Université, Collège de France, 24 rue Lhomond, 75005 Paris, France

###### Abstract

Three-photon (3P) microscopy is getting traction due to its superior
performance in deep tissues. Yet, aberrations and light scattering still pose
one of the main limitations in the attainable depth ranges for high-resolution
imaging. Here, we show scattering correcting wavefront shaping with a simple
continuous optimization algorithm, guided by the integrated 3P fluorescence
signal. We demonstrate focusing and imaging behind scattering layers and
investigate convergence trajectories for different sample geometries and
feedback non-linearities. Furthermore, we show imaging through a mouse skull
and demonstrate a novel fast phase estimation scheme that substantially
increases the speed at which the optimal correction can be found.

Light passing through biological tissue experiences an environment of
heterogeneous refractive index, inevitably leading to aberrations and
scattering Ntziachristos (2010). These effects limit optical microscopy to
superficial tissue layers. Adaptive optics (AO), originally developed for
astronomy, allows to correct for common lower order aberrations, substantially
improving imaging quality at depth Hampson  _et al._ (2021). In addition,
point-scanning techniques that provide intrinsic optical sectioning, like
confocal or two-photon (2P) microscopy, are likewise pushing the depth limit
further into the tissue Wilson (2011); Helmchen and Denk (2005).

Still, light scattering poses one of the ultimate limitations for imaging deep
inside living tissues. It is, however, a coherent and deterministic process
that can be controlled, even in regimes where no ballistic light is present
Vellekoop and Mosk (2007). To form a focus inside such strongly scattering
environments non-invasively, one needs a metric providing feedback on the
quality of the focus within the sample Horstmeyer  _et al._ (2015). In multi-
photon fluorescence imaging, such a metric is conveniently given by the total
fluorescence signal, due to the non-linear nature of the excitation process: a
tight focus on a single target generates more signal than distributing the
same power over multiple targets. For 2P fluorescence (2PF), this feedback has
been shown to facilitate focusing and imaging for a broad set of wavefront
optimization algorithms Katz  _et al._ (2014); Papadopoulos  _et al._ (2017);
May  _et al._ (2021); Blochet  _et al._ (2021).

Recently, three-photon (3P) microscopy got adopted for biological imaging,
primarily due to its superior performance at depth Horton  _et al._ (2013);
Wang  _et al._ (2018); Weisenburger  _et al._ (2019). Both, the higher non-
linearity of the excitation process and the shift towards longer excitation
wavelength extends the achieved depth ranges. Nevertheless, AO still proves
highly useful to improve the 3PF signal when imaging at depth Rodríguez  _et
al._ (2021); Streich  _et al._ (2021); Sinefeld  _et al._ (2022); Qin  _et
al._ (2022). This is because the higher non-linearity also leads to a higher
susceptibility to signal loss due to aberrations. Going beyond the correction
of low-order aberrations is therefore hard, due the vanishing initial signal
obtained from a strongly speckled point-spread function (PSF). Even though,
for particular optimization schemes, focusing and imaging has been
demonstrated in such a context Berlage  _et al._ (2021), a general
understanding of the convergence criteria and optimal strategies for wavefront
shaping with 3PF feedback are still missing.

Here, we demonstrate blind focusing behind a scattering layer through a simple
continuous wavefront optimization scheme Vellekoop and Mosk (2008), informed
only by the total 3PF signal. Scanning the obtained focus within the range of
the memory effect Feng  _et al._ (1988), we show signal improvements of up to
two orders of magnitude and diffraction-limited imaging in situations where
the unshaped PSF only resolves large-scale features. Our approach reliably
achieves focusing, also in volumetrically labeled three-dimensional (3D)
samples. For volumetrically labeled structures, it was recently argued that
the enhanced non-linearity of 3P excitation is a prerequisite for convergence
and that 2PF feedback could not achieve convergence to a focus in 3D samples
Katz  _et al._ (2014); Berlage  _et al._ (2021). We therefore compare the
convergence of continuous wavefront optimizations guided by 3PF and 2PF
feedback for different sample geometries and show that 2PF feedback does
indeed suffice to focus in volumetric targets. We further demonstrate focusing
and imaging through a thinned mouse skull bone, a scattering environment
particularly relevant for neuroscience applications. Finally, we implement a
novel two-point phase estimation procedure, achieving a more than twofold
improvement in speed for continuous wavefront optimization schemes,
highlighting the practicality of this simple approach.

Figure 1: Continuous wavefront optimization guided by 3PF feedback. (a) Scheme
of the microscope setup. In the top right, the fluorescence signal recorded on
the PMT when phase stepping a single mode is shown. Examples of circular
Hadamard patterns are depicted in the bottom right. The central zoom-in shows
the bead sample above the 2D scattering surface. (b) Evolution of the
fluorescence signal enhancement during optimization. The gray and white
backgrounds indicate a full set of 256 modes. (c) Final correction pattern
applied to the SLM. (d) PSF of the unshaped excitation beam, measured in
transmission. (e) Image of the sample obtained from scanning the unshaped PSF
shown in (d) over the sample. (f) PSF after wavefront shaping. (g) The
corresponding scanned image.

We use a custom-built microscope with a liquid crystal based spatial light
modulator (SLM) conjugated to the back focal plane of the illumination
objective (Fig. 1(a), for details see SM Sect. I). A second objective aligned
in transmission collects the emitted fluorescence which is measured on a
photo-multiplier tube (PMT). A CCD camera in the transmission path further
enables monitoring of the excitation light pattern at the sample plane during
the wavefront shaping procedure. As samples, we use fluorescent beads which
are $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in diameter. The beads are drop-
cast on a glass coverslip whose bottom surface, facing the illumination
objective, was sand-blasted to introduce a single strongly scattering layer
(see SM Sect. II). This ensures a large speckle envelope and the absence of
any ballistic light.

To optimize the incident wavefront we use a continuous optimization algorithm
Vellekoop and Mosk (2008). In each optimization step, half of the active
pixels on the SLM are phase-stepped between 0 and 2$\pi$ while the other half
is kept constant and acts as a reference. From the observed 3PF signal, we
extract the phase at which the modulated pixels are optimal and add it to the
corresponding pixels of the current optimal pattern. The resulting pattern
forms the initial state of the next optimization step and is thereby
continuously improved. The mode patterns, determining the active pixels in
each optimization step, are chosen to be circular Hadamard patterns (Fig. 1(a)
bottom right, for details see SM Sect. III). While more involved procedures
have been used for focusing through wavefront shaping with non-linear feedback
Papadopoulos  _et al._ (2017); May  _et al._ (2021), this algorithm represents
one of the most simple ones.

A typical evolution of the 3PF signal under optimization is shown in Fig.
1(b). Here, the active area on the SLM was decomposed into 256 modes which
were optimized three times, consecutively. We see that already after the
second round the fluorescence signal almost plateaued to its maximal value, a
number of measurements comparable to other approaches Berlage  _et al._
(2021). The final SLM pattern (Fig. 1(c)) reaches a 60-fold signal
enhancement, realizing a tight focus of the excitation light at the target
plane (Fig. 1(f)). Note that, while the maximal intensity of the excitation
light increases by about a factor of 10, the fluorescence signal enhancement
stays way below the cube of this value. This is because the target sample is
extended such that the initial speckled PSF excites many beads at once and the
redistribution of power into a single focus reduces the signal collected from
neighboring beads.

Applying additional linear phase ramps to the SLM allows us to scan the beam
over the sample without dedicated scanning mirrors, albeit at a lower speed.
Within the memory effect range of the scattering geometry, the optimized
correction pattern stays valid and we can form an image of the hidden sample.
Additionally, laterally translating the correction pattern on the SLM
proportionally to the phase gradients applied for each position extends the
validity range of the correction further Osnabrugge  _et al._ (2017);
Papadopoulos  _et al._ (2020) (see SM Sect. IV). The resulting image with the
beads well resolved can be seen in Fig. 1(g). In contrast, in an image formed
by scanning the unshaped PSF over the sample we could only resolve global
shapes (Fig. 1(d) and (e)). This is different to many multi-photon AO
implementations where often a substantial signal enhancement is achieved while
improvements in resolution remain small.

Due to the cubic non-linearity of the feedback signal and the multitude of
sources illuminated by the initial speckle pattern it is not a given that the
optimization converges to a single focus. Yet, we observe reliable focusing on
similar time scales for many different speckles and sample geometries, even
though the exact position of the focus cannot be predicted a priori. The
question how the convergence to a focus depends on the non-linearity of the
feedback process and the geometry of the sample is a point of interest in
recent studies on this topic. It was argued that for volumetrically labeled
samples 2PF feedback does not suffice to form a single focus Katz  _et al._
(2014); Berlage  _et al._ (2021). The argument is based on the observation
that under 2P excitation the signal collected from a Gaussian excitation beam
inside a volumetric dye pool does not depend on the NA of the beam. Only for a
third order non-linearity this dependence appears. Nevertheless, in the
context of 2P AO, convergence was observed for aberrated beams in volumetric
samples Sinefeld  _et al._ (2015); Tang  _et al._ (2012). For a 2D geometry,
it was also predicted that for certain optimization schemes there exists a
maximal number of targets such that wavefront shaping under non-linear
feedback reliably converges to a focus Osnabrugge  _et al._ (2019).

Figure 2:  Convergence study. (a) Schematic of the sample geometries. 1D: an
isolated bead $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in diameter, 2D: a
homogeneous layer of the same beads, 3D: a large bead
$15\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in diameter. (b) The upper panel
shows a comparison of signal enhancement during wavefront shaping for
different geometries under 2P feedback. The gray and white areas indicate
single repetitions of the set of 256 modes used. The lower panel shows the
same data in log scale. (c) The same for 3P feedback. All curves are averaged
over at least four realizations of the scattering layer.

We therefore investigate the convergence properties of the wavefront shaping
algorithm for different sample geometries under 2P and 3P excitation. For
that, we prepare samples with an effective 3D, 2D, and 1D distribution of
fluorophores (Fig. 2(a)). The 1D samples consist of a single bead while 2D
samples are created from bead monolayers which allows for a simple preparation
with a precise thickness. For the 3D samples we use large beads with a
diameter of $15\text{\,}\mathrm{\SIUnitSymbolMicro m}$. Being about the size
of a neuron soma, such finite volumetric targets are much more relevant to
biomedical applications than dye pools. Also, when interested in geometric
effects, they assure that convergence is not hindered by the total number of
targets Osnabrugge  _et al._ (2019).

Figure 2(b) and (c) show the signal enhancement under optimization for 2P and
3P excitation, respectively. We observe that for 2P excitation, the difference
in enhancements between 1D and extended targets is much larger than in the
case of 3P excitation. We attribute this to the geometric effects discussed
above. Interestingly, comparing the curves for the different geometries in log
scale, we see no large differences in convergence time within one excitation
regime. In all cases, the final signal level is reached after about three
optimization runs through the full mode basis.

Figure 3:  Focusing through a thinned mouse skull bone. (a) Initial speckled
PSF before wavefront shaping. (b) Focus obtained after optimizing 256 modes,
each two times. (c) Image obtained from scanning the PSF shown in (a) over the
bead sample. (d) Corresponding image formed by scanning the shaped PSF.

We further observe that 2P feedback does indeed lead to a single focus within
the large bead volume. Even though the signal enhancement is comparatively
small, the algorithm converges reliably. This points to the conclusion that
there is no qualitative difference for 2P feedback in volumetric samples. Yet,
as the increase per mode is small, for some practical situations, noise levels
can render 2P feedback ineffective. While it could be argued that the finite
volume used here aids convergence, we also observe focusing in dye pool
samples (see SM Sect. V). Note that, even though the small 2P signal increase
for volumetric samples renders wavefront shaping ineffective in weakly
scattering regimes, in the case of a fully speckled PSF with no strong
ballistic component the situation is different. There, correcting the
wavefront enables imaging, even if the initial uncorrected signal is of
approximately the same magnitude

To show the practical applicability of this wavefront shaping technique we use
it to focus through a relevant biological obstacle: a mouse skull bone. The
skull was thinned to about $200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
fluorescent beads were placed behind it. Figure 3(a) and (b) show a PSF
measured in transmission before and after shaping, demonstrating focusing
after twice optimizing 256 modes. The images obtained with these PSFs show a
similar improvement as for the single scattering layer, albeit with a slightly
smaller total enhancement and a more limited memory effect range (Fig. 3(c)
and (d)). This is expected as the skull bone is 3D scattering obstacle,
compared to the 2D ground glass layer.

Figure 4:  Fast phase estimation. Evolution of the fluorescence signal
enhancement for an optimization using 6 phase steps for the initial 128 modes
and 2 phase steps thereafter. The two insets show examples of phase estimation
(blue dots) for the 6-point (top left) and 2-point measurements (bottom
right). The orange lines show the expected non-linear fluorescence signal with
the phase extracted from the corresponding methods.

The downside of a modulation scheme where in each iteration half the SLM
pixels are phase stepped while the other half acts as a reference is the
pronounced non-linearity of the feedback signal. Even though it leads to the
maximal signal contrast, this configuration requires more than 4 measurements
per mode, especially for 3rd order non-linearities or above (see SM Sect. VI).
For the measurements presented so far we therefore used 6 measurement points
to determine the optimal phase for each mode. However, when monitoring the
estimated phases during the optimization one realizes that they quickly tend
to zero after the first few iterations. This is because the applied wavefront
starts to converge to the optimal one, rendering each further improvement
incremental. Knowing this we can devise an optimal measurement strategy that
does not require a full 6 measurements for each iteration (see SM Sect. VII).
Choosing two measurement points that coincide with the maximal gradient of the
non-linear fluorescence signal allows us to estimate the optimal phase using
far less measurements. This is reminiscent of the standard 3-point
measurements scheme used in AO Debarre  _et al._ (2007); Sinefeld  _et al._
(2015), only that here we know the functional shape of the signal and can
precisely determine the optimal stroke with just two data points. Figure 4
shows the application of this technique, demonstrating a speedup of more than
a factor of two for the full optimization procedure. Note that, at the point
in the optimization where the fast 2-point phase estimation starts to be
applied, the signal is only 3% of its later maximal value.

Recently, several novel wavefront shaping strategies aimed at focusing inside
scattering environments have been developed Papadopoulos  _et al._ (2017); May
_et al._ (2021). Here, we demonstrate that for 3PF samples a simple continuous
wavefront shaping procedure suffices to blindly focus behind a scattering
layer, even in the absence of any initial ballistic component. It reliably
converges and can even be run at two measurements per mode, after some initial
pre-optimization. We show imaging in situations where initially only large
structures can be resolved and further demonstrate its application with
biologically relevant scattering media. Investigating the convergence in
continuous optimization for different target geometries using both 2PF and 3PF
feedback we further demonstrate the advantage a higher non-linearity gives in
the signal increase obtained for extended sources. We also show that, contrary
to the intuition obtained from a Gaussian beam model, 2PF feedback does allow
focusing inside homogeneous volumetric targets, albeit with low signal
enhancements. This highlights that wavefront optimization based on non-linear
feedback is far from trivial and we hope the results presented here inspire
further numerical and experimental studies.

Acknowledgments:
This project was funding by the European Research Council under the grant
agreement No. 724473 (SMARTIES), the European Union’s Horizon 2020 research
and innovation program under the FET-Open grant No. 863203 (Dynamic) and the
program ”Investissements d’Avenir” launched by the French Government and
implemented by ANR with the references ANR–10–LABX–54 MEMOLIFE and
ANR–10–IDEX–0001–02 PSL* Université Paris. B.R. was supported by the European
Union’s Marie Skłodowska-Curie fellowship, grant agreement No. 888707
(DEEP3P).

The authors thank Walther Akemann for providing the mouse skull; Julien
Guilbert and Gerwin Osnabrugge for useful discussion and and Fei Xia for the
careful reading of the manuscript.

## References

  * Ntziachristos (2010) V. Ntziachristos, Nature Methods 7, 603 (2010).
  * Hampson  _et al._ (2021) K. M. Hampson, R. Turcotte, D. T. Miller, K. Kurokawa, J. R. Males, N. Ji, and M. J. Booth, Nature Reviews Methods Primers 1, 68 (2021).
  * Wilson (2011) T. Wilson, Journal of Microscopy 244, 113 (2011).
  * Helmchen and Denk (2005) F. Helmchen and W. Denk, Nature Methods 2, 932 (2005).
  * Vellekoop and Mosk (2007) I. M. Vellekoop and A. P. Mosk, Optics Letters 32, 2309 (2007).
  * Horstmeyer  _et al._ (2015) R. Horstmeyer, H. Ruan, and C. Yang, Nature Photonics 9, 563 (2015).
  * Katz  _et al._ (2014) O. Katz, E. Small, Y. Guan, and Y. Silberberg, Optica 1, 170 (2014).
  * Papadopoulos  _et al._ (2017) I. N. Papadopoulos, J.-S. Jouhanneau, J. F. A. Poulet, and B. Judkewitz, Nature Photonics 11, 116 (2017).
  * May  _et al._ (2021) M. A. May, N. Barré, K. K. Kummer, M. Kress, M. Ritsch-Marte, and A. Jesacher, Nature Communications 12 (2021), 10.1038/s41467-021-24666-9.
  * Blochet  _et al._ (2021) B. Blochet, W. Akemann, S. Gigan, and L. Bourdieu, bioRxiv (2021), 10.1101/2021.09.06.459064.
  * Horton  _et al._ (2013) N. G. Horton, K. Wang, D. Kobat, C. G. Clark, F. W. Wise, C. B. Schaffer, and C. Xu, Nature Photonics 7, 205 (2013).
  * Wang  _et al._ (2018) T. Wang, D. G. Ouzounov, C. Wu, N. G. Horton, B. Zhang, C.-H. Wu, Y. Zhang, M. J. Schnitzer, and C. Xu, Nature Methods 15, 789 (2018).
  * Weisenburger  _et al._ (2019) S. Weisenburger, F. Tejera, J. Demas, B. Chen, J. Manley, F. T. Sparks, F. Martínez Traub, T. Daigle, H. Zeng, A. Losonczy, and A. Vaziri, Cell , 1 (2019).
  * Rodríguez  _et al._ (2021) C. Rodríguez, A. Chen, J. A. Rivera, M. A. Mohr, Y. Liang, R. G. Natan, W. Sun, D. E. Milkie, T. G. Bifano, X. Chen, and N. Ji, Nature Methods 18, 1259 (2021).
  * Streich  _et al._ (2021) L. Streich, J. C. Boffi, L. Wang, K. Alhalaseh, M. Barbieri, R. Rehm, S. Deivasigamani, C. T. Gross, A. Agarwal, and R. Prevedel, Nature Methods 18, 1253 (2021).
  * Sinefeld  _et al._ (2022) D. Sinefeld, F. Xia, M. Wang, T. Wang, C. Wu, X. Yang, H. P. Paudel, D. G. Ouzounov, T. G. Bifano, and C. Xu, Frontiers in Neuroscience 16, 1 (2022).
  * Qin  _et al._ (2022) Z. Qin, Z. She, C. Chen, W. Wu, J. K. Y. Lau, N. Y. Ip, and J. Y. Qu, Nature Biotechnology (2022), 10.1038/s41587-022-01343-w.
  * Berlage  _et al._ (2021) C. Berlage, M. Tantirigama, M. Babot, D. D. Battista, C. Whitmire, I. N. Papadopoulos, J. F. A. Poulet, M. Larkum, and B. Judkewitz, bioRxiv 8, 1 (2021).
  * Vellekoop and Mosk (2008) I. Vellekoop and A. Mosk, Optics Communications 281, 3071 (2008).
  * Feng  _et al._ (1988) S. Feng, C. Kane, P. A. Lee, and A. D. Stone, Physical Review Letters 61, 834 (1988).
  * Osnabrugge  _et al._ (2017) G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, and I. M. Vellekoop, Optica 4, 886 (2017), arXiv:1705.01373 .
  * Papadopoulos  _et al._ (2020) I. N. Papadopoulos, J. S. Jouhanneau, N. Takahashi, D. Kaplan, M. Larkum, J. Poulet, and B. Judkewitz, Light: Science and Applications 9 (2020), 10.1038/s41377-020-00348-x.
  * Sinefeld  _et al._ (2015) D. Sinefeld, H. P. Paudel, D. G. Ouzounov, T. G. Bifano, and C. Xu, Optics Express 23, 31472 (2015).
  * Tang  _et al._ (2012) J. Tang, R. N. Germain, and M. Cui, Proceedings of the National Academy of Sciences 109, 8434 (2012).
  * Osnabrugge  _et al._ (2019) G. Osnabrugge, L. V. Amitonova, and I. M. Vellekoop, Optics Express 27, 11673 (2019).
  * Debarre  _et al._ (2007) D. Debarre, M. J. Booth, and T. Wilson, Optics Express 15, 8176 (2007).

Supplementary Materials

##  I Experimental setup

The experimental setup is sketched in Fig. 1(a) of the main text. The source
used is an OPO (Coherent Mira OPO-X) pumped by a pulsed Ti-sapphire laser
(Coherent Chameleon Ultra II) with $140\text{\,}\mathrm{fs}$ pulse duration
and repetition rate of $80\text{\,}\mathrm{MHz}$. To modulate the wavefront,
we use a liquid crystal phase-only SLM (Meadowlarks HSP512L) which is imaged
on the back focal plane of a water immersion objective (Zeiss W Plan-
Apochromat 40x/1.0 DIC M27). After illuminating the sample, a second
microscope objective (Zeiss EC Plan-NEOFLUAR 40x/1.3) images the excitation
light pattern onto a CCD camera (Basler acA1300-30um) and collects the
fluorescence light which is diverted towards a PMT (Hamamatsu H7422P-40) by a
dichrocic mirror (Semrock Di03-R785-t1-25x36). Additional filters in front of
the PMT (Thorlabs FESH0600, FESH0700) assure that only fluorescent light is
detected.

For 3P excitation, the OPO is tuned to $1050\text{\,}\mathrm{nm}$. As targets,
we use fluorescent beads with a single photon absorption peak in the blue that
are $2.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in diameter (Thermo-Fisher
Fluoro-max B0200 - 365, 388, $412\text{\,}\mathrm{nm}$ / 445, 445,
$473\text{\,}\mathrm{nm}$). The large beads used as 3D targets for the
measurements presented in Fig. 2 of the main text are
$15\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in diameter and have a similar
single photon absorption peak (Invitrogen F8837 - $365\text{\,}\mathrm{nm}$ /
$415\text{\,}\mathrm{nm}$). For the 2P excitation measurements presented in
Fig. 2 of the main text, we use the same beads excited directly by the light
from the pump laser at $800\text{\,}\mathrm{nm}$.

Typically, for 3P excitation, low repetition rate lasers are used, as for a
repetition rate of $80\text{\,}\mathrm{MHz}$, the 3P the absorption is quite
inefficient. We therefore use initial laser powers of about
$180\text{\,}\mathrm{mW}$ at the sample plane and detect the fluorescence in
transmission with a high NA oil immersion objective. Note that, the choice
against detecting fluorescence in an epi-geometry is solely based on the need
to increase signal levels and is no requirement of the optimization strategy
nor the imaging technique used. In addition, we underfill the illumination
objective back aperture to prevent clipping and the associated power loss. The
Gaussian beam entering the aperture has a $1/e^{2}$ width of
$6.3\text{\,}\mathrm{mm}$ for the 3P excitation at $1050\text{\,}\mathrm{nm}$
and $8\text{\,}\mathrm{mm}$ for the 2P excitation at
$800\text{\,}\mathrm{nm}$. The back aperture of the illumination objective is
about $9\text{\,}\mathrm{mm}$ wide. On the SLM, we modulate a circular region
(see Sect. III) that spans about 80% of the $1/e^{2}$ beam width, for both 3P
and 2P excitation.

To exclude the correction of aberrations inherent to the optical setup, we
measure the optimal system correction pattern and use it as the initial point
for the wavefront optimization process. For that, we use samples with a single
fluorescent bead of diameter $0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$
(Invitrogen F8812 and F8813). These bead samples are prepared in the same way
as the samples used for the measurements presented in the main text (see Sect.
II) but without a scattering layer. Starting from a flat SLM pattern and
optimizing the total non-linear fluorescent signal when focusing on one of
these beads we obtain the system correction pattern. For higher signal
strength, we use a 2P excitation for both the $1050\text{\,}\mathrm{nm}$ and
the $800\text{\,}\mathrm{nm}$ light. The final pattern serves as the initial
wavefront for the optimization procedures presented in the main text. Note
that, the signal increase due to the system correction is maximally 25% and
therefore negligible compared to the one achieved in scattering correction.

In our setup, the speed of the wavefront optimization is limited by the
integration time applied to the PMT signal, which is set to
$40\text{\,}\mathrm{ms}$. Again, this is necessary due to the low signal
levels. With the $10\text{\,}\mathrm{ms}$ it takes to load a new SLM pattern
and some processing time and additional delays that assure smooth operation,
this results in $70\text{\,}\mathrm{ms}$ for a single data point and
$420\text{\,}\mathrm{ms}$ to perform a 6-point phase stepping measurement. In
the case of a 2-point phase estimation, as presented in Fig. 4 of the main
text, the time to measure one mode reduces to $140\text{\,}\mathrm{ms}$. Note
that, using a source more adapted to 3P excitation would reduce this time
substantially.

The signal increase during wavefront shaping can be large, up to 280-fold in
the case of a single bead under 3P excitation, as demonstrated in Fig. 2(c) of
the main text. To prevent bleaching and saturation of the PMT we therefore
reduce the laser power during the optimization. This is automated with a
continuously variable neutral density filter (Thorlabs NDC-50C-4) rotated by a
computer controlled stepper motor. The power is reduced as soon as the
threshold of a 7 to 10-fold signal increase is reached.

##  II Sample preparation

A scattering surface is introduced to a #1.5 cover slip by blasting its bottom
surface with 220 grit sand. To reduce the scattering, the resulting rough
surface is covered with Norland NOA65 glue which is immediately wiped off with
a lens tissue. After that, the glue remaining on the surface is cured under UV
light. The beads are drop-casted on the top surface of the coverslip, about
$170\text{\,}\mathrm{\SIUnitSymbolMicro m}$ above the scattering layer. They
are covered by the same NOA65 glue to create a near-index matched environment.
Finally, a second #1.5 cover slip is placed on top of the beads and the glue
is cured. This second curing process is kept short not to bleach the beads.
Near-index matching the polystyrene beads ensures that the excitation light
recorded in transmission is not distorted by the beads.

For the measurement presented in Fig. 3 of the main text, a piece of the
parietal bone of a male 7.5 weeks old C57BL/6 mouse was used. The bone was
thinned down to $\sim$$200\text{\,}\mathrm{\SIUnitSymbolMicro m}$, immersed in
phosphate-buffered saline and sandwiched between two #1.5 cover slips
separated by a $240\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick spacer. The
beads are drop-cast on the top surface of the upper coverslip and covered by
UV-cured glue in the same procedure as described above.

Experimental procedures were conducted in accordance with the institutional
guidelines and in compliance with French and European laws and policies. They
were approved by the ‘Charles Darwin’ local institutional ethical committee
registered at the French National Committee of Ethical Reflection on Animal
Experimentation under the number 05 (authorization number: APAFIS 26667)

##  III Circular Hadamard modes

The rows of the Hadamard matrix, reshaped into 2D patterns, are a common
choice of basis for wavefront shaping applications beyond Zernike polynomials.
They are orthogonal and each mode divides the pixels in the modulated area
into two equally large sets. The latter property assures a maximal
interference signal when modulating one of the sets against the other (see
Sect. VI) while the former provides a faster convergence in an iterative
optimization scheme compared to, e.g., a random basis.

In microscopy, however, we often want to control the light field at the
backfocal plane of an objective or at some other circular symmetric optical
component. Here, the square shaped Hadamard modes are sometimes unpractical as
they can be clipped by back apertures or lead to an asymmetric modulation of
the light. The square shape though is arbitrarily chosen and not inherent to
the originally 1D rows of the Hadamard matrix. We can arrange the rows in any
kind of 2D pattern while still preserving the patterns’ orthogonality and
equal-division of pixels. A circular arrangement for example, as shown in Fig.
S1(a), can be optimal to when working with circular optics and apertures. The
resulting mode patterns are slightly wider than their square equivalents, as
shown for two examples in Fig. S1(b).

Figure S1:  Circular Hadamard modes. (a) Circular pixel layout for
$N_{\mathrm{mode}}=256$. The rows of the Hadamard matrix are arranged
following the white arrow. (b) Four examples of square Hadamard modes and
their circular counterparts (pairs left and right, respectiveley).

Note that, while for $N_{\mathrm{mode}}=256$ a perfect circle can be arranged,
for other basis set sizes a few pixels have to be added or removed
asymmetrically to assure a number of pixels equal to a power of two, as
required by the Hadamard matrix. These few pixels, however, generally do not
affect the circular shape much.

##  IV Shift scanning

To scan the optimized focus over the sample, forming an image, we apply phase
ramps to the SLM. This is a very slow imaging technique which could be easily
sped up by introducing galvanometric mirrors to the beam path after the SLM.
However, performing the scan with the SLM provides the advantage of being able
to correct the applied shaping pattern while scanning. For example, one can
partially correct the fact that while scanning, the scattering regions the
beam traverses shift with respect to beam, misaligning the correction pattern
and the scatterers Osnabrugge  _et al._ (2017); Papadopoulos  _et al._ (2020).
Simply shifting the correction pattern on the SLM as a function of the applied
phase gradient leads to a substantially enlarged memory effect range, as shown
in Fig. S2. Here, the optimal amount of shift is found iteratively, by tuning
the parameter while observing the range at which the wavefront correction is
valid. For the images presented in the main text a shift per phase gradient of
$0.032\text{\,}\mathrm{m}\mathrm{m}^{2}\mathrm{/}\mathrm{rad}$ was used. At
the edge of the image this corresponds to a 19 pixel shift, which is about 10%
of the total pattern length. With the thin 2D scattering layers, we observed
an extension of the memory effect range of excitation intensity by a factor of
2.7. For volumetric scatterers, like the mouse skull bone of Fig. 3 of the the
main text, the gain in imaging range is smaller.

Figure S2:  Shift-enhanced memory effect range. (a) 2PF image obtained from
scanning a speckled PSF over a bead target. (b) Image obtained after wavefront
shaping, applying phase gradients to the SLM to move the focus over the
sample. (c) Same image when, in addition to the phase gradients, the optimal
correction mask is shifted in the direction of the gradient.

##  V Dye pool measurement

To assure that the convergence to a focus with 2PF feedback observed in the
measurements presented in Fig. 2 of the main text is not due to the finite
size of the fluorescent target we also performed the same measurement on a dye
pool sample with a thickness of $120\text{\,}\mathrm{\SIUnitSymbolMicro m}$.
The results are presented in Fig. S3. The initial speckle converges towards a
clear focus although the increase in total fluorescence is only 70%. Also, the
convergence is slower than observed in the measurements of the main text, even
after sampling through all modes four times the wavefront does not seem to
have fully converged.

Figure S3:  Focusing guided by 2PF feedback inside a dye pool. (a) PSF of the
unshaped excitation beam, measured in transmission. (b) PSF after wavefront
shaping. (c) Evolution of the fluorescence signal enhancement during
optimization. The gray and white backgrounds indicate a full set of 256 modes.
The enhancement is averaged over three scattering realizations.

Note that, for these measurements we needed to reduce the size of the initial
speckle envelope to observe convergence. We believe this is related to the
limits encountered when exciting a large number of targets, discussed in
Osnabrugge  _et al._ (2019), although 3D targets are not discussed there.

To reduce the scattering and thereby shrink the envelope we use a glue at a
lower refractive index (Norland NOA85) and fully cover the sand blasted cover
glass surface without wiping it off. We then add an additional coverslip on
the bottom and cure the glue. This reduces the refractive index difference at
the rough interface, limiting the scattering to smaller angles.

##  VI Phase estimation

Tuning the phase of a given mode during the optimization procedure leads to a
total $n$-photon fluorescent signal

$F(\delta)\propto\sum_{i}\sigma_{i,n}\big{[}\lvert
E^{i}_{\mathrm{ref}}\rvert^{2}+\lvert E^{i}_{\mathrm{mode}}\rvert^{2}+2\lvert
E^{i}_{\mathrm{ref}}\rvert\lvert
E^{i}_{\mathrm{mode}}\rvert\mathrm{cos}(\varphi_{i}+\delta)\big{]}^{n}.$ (1)

Here, $E^{i}_{\mathrm{mode}}$ is the field of the tuned mode at the point of
source $i$ while $E^{i}_{\mathrm{ref}}$ is the reference field at that same
point. The phase between these two fields is given by $\varphi_{i}$ and
$\delta$ denotes the applied phase stepping. The signal is detected on a
bucket detector such that the sum runs over all sources $i$ that are
illuminated, with $\sigma_{i,n}$ being the $n$-photon cross section of each
source.

Typically, in continuous wavefront optimization schemes, the optimal phase of
a single mode is determined by stepping $\delta$ from 0 and $2\pi$ at equal
spacing. Assuming an offset cosine modulated signal and enough measurement
points, the optimal phase is then given by

$\varphi_{\mathrm{opt}}=\mathrm{angle}\Big{[}\sum_{\delta}F(\delta)e^{i\delta}\Big{]}.$
(2)

However, the cosine shape of $F(\delta)$ is only assured for linear feedback
mechanisms. For $n>1$, the functional form of the total fluorescent signal
depends on the target and the amplitude balance between reference and mode
fields. Expanding the $n$-th power in eq. 1, we can see that only if one of
the two fields is much larger than the other a term linear in
$\mathrm{cos}(\varphi_{i}+\delta)$ dominates the $\delta$ dependent terms.
Therefore, it is a viable strategy to choose $\lvert
E^{i}_{\mathrm{ref}}\rvert\ll\lvert E^{i}_{\mathrm{mode}}\rvert$ or vice versa
in order to assure a cosine shaped modulation signal. The downside is a small
modulation amplitude compared to the background fluorescence.

For $\lvert E^{i}_{\mathrm{ref}}\rvert=\lvert E^{i}_{\mathrm{mode}}\rvert$ the
modulation amplitude is the largest, but for a single target
$F(\delta)\propto\big{[}1+\mathrm{cos}(\varphi_{\mathrm{opt}}+\delta)\big{]}^{n}$
and for multiple fluorescent sources the superposition of many non-linear
modulations leads to strongly target dependent functional forms. Yet, eq. 2
still provides a good estimate of $\varphi_{\mathrm{opt}}$ in this case and a
continuous wavefront optimization based on it converges quickly. A critical
point, however, is the number of phase steps acquired for each mode. While
typically a 4-point phase stepping is used, for non-linear fluorescence
signals this can be insufficient, as shown in Fig. S4(a) for a 3PF signal from
a single target. For certain optimal phases, the measured points are sub-
optimally placed in regions of flat or slowly varying signals leading to
increased errors in the estimation of $\varphi_{\mathrm{opt}}$. This can be
counteracted by measuring $F(\delta)$ at more points, as shown in Fig. S4(b),
with the downside of a reduced optimization speed.

Figure S4: Phase estimation for continuous wavefront shaping. (a) Phase
estimation with 4 equally spaced measurement points (blue dots) for a
maximally non-linear feedback signal originating from a single 3PF source
(orange line). (b) The same with 6 equally spaced measurement points. (c)
Simulated distribution of $\varphi_{\mathrm{opt}}$ for a continuous wavefront
optimization. The scattering is modeled by a complex Gaussian i.i.d.
transmission matrix with 256 input modes. This set of modes is optimized two
times with the dashed vertical line indicating the two sets. The 3PF target is
extended and covers 27% of the output pixels. (d) Amplitude ratio of the
currently optimized mode and the reference field at the focal point for the
same simulation. (e) Phase estimation with 2 measurement points at
$\tilde{\delta}_{\pm}$. The dash-dotted lines with slopes $\pm\;\alpha/2$
indicate the first order phase estimation. The orange circles show the large
signal changes even for $\varphi_{\mathrm{opt}}\simeq 0$. For (a), (b) and (e)
$\varphi_{\mathrm{opt}}=0.05\;\pi$.

For a continuous wavefront optimization procedure using Hadamard modes where
in each step half of the SLM pixels are optimized and the other half acts as
the reference, the two interfering fields are on average balanced leading to
non-cosine modulations of the total fluorescence signal. Therefore, for the
data presented in Fig. 1, 2 and 3 of the main text, a 6-point optimization
scheme is used.

Another interesting point to note is that for a continuous optimization
strategy where always half of the pixels are modulated, the reference field
changes in each iteration. This leads to a slower conversion, as pixels are
constantly optimized relative to a different reference. Even in case of a
linear feedback signal, such an optimization does not fully converge after
sampling through all modes once. The requirement to optimize each mode 2-3
times, observed in Fig. 1(b), Fig. 2 and Fig. 4 of the main text, is therefore
not solely a result of the non-linear feedback signal.

##  VII 2-point phase estimation

For continuous wavefront optimization schemes, the current wavefront gets
closer and closer to one that focuses on a single target as the optimization
continues. Also, the partition in the modulated part of the wavefront and the
reference part changes in each iteration. This means, as the optimization
continuous, both the current mode and the reference field tend to be already
close to their optimal orientation even before the phase stepping procedure
starts. We can observe this when looking at the distribution of the extracted
optimal phases or when comparing the amplitudes of the current mode and the
reference at the brightest spot during the optimization (see Fig. S4(c) and
Fig. S4(d), respectively, for simulated distributions). Together with the
observation that even for complicated and large targets the procedure quickly
converges to a single focus which is then slowly optimization, this leads to
predictable fluorescence signal after a few iterations. The balance between
$\lvert E_{\mathrm{ref}}\rvert$ and $\lvert E_{\mathrm{mode}}\rvert$ forces
$F(\delta)\propto\big{[}1+\mathrm{cos}(\varphi_{\mathrm{opt}}+\delta)\big{]}^{n}$,
as discussed above, and $\varphi_{\mathrm{opt}}$ is always close to zero.

This knowledge can be used to selectively measure only at those phase points
that provide the most information and thereby speed up the optimization. From
Fig. S4(b) for example, we see that the two points in the flanks of the peak
will be most crucial in determining $\varphi_{\mathrm{opt}}$. Ideally, we want
to measure at the maximal signal gradient, which occurs at

$\bigg{[}\frac{1+\mathrm{cos}(\tilde{\delta})}{2}\bigg{]}^{n}=\frac{1}{2},$
(3)

determining the optimal measurement points $\tilde{\delta}_{\pm}$ to be given
by

$\tilde{\delta}_{\pm}=\pm\;\mathrm{arccos}(2^{\frac{n-1}{n}}-1).$ (4)

Expanding $F(\delta)$ around $\varphi_{\mathrm{opt}}=0$ up to first order
further provides us with a compact linear expression to estimate
$\varphi_{\mathrm{opt}}$

$\varphi_{\mathrm{opt}}\simeq\frac{1}{\alpha}\;\frac{F(\tilde{\delta}_{+})-F(\tilde{\delta}_{-})}{F(\tilde{\delta}_{+})+F(\tilde{\delta}_{-})}$
(5)

with

$\alpha=n\;\mathrm{sin}(\tilde{\delta}_{+})\;\bigg{[}\frac{1+\mathrm{cos}(\tilde{\delta}_{+})}{2}\bigg{]}^{n-1}.$
(6)

A graphical illustration of the 2-point phase estimation is shown in Fig.
S4(e). After starting the optimization procedure with regular equally spaced
phase stepping, once $\varphi_{\mathrm{opt}}$ starts to converge to small
angles the 2-point measurement scheme can be employed. In the measurements
presented in Fig. 4 of the main text, the 2-point scheme was started after
$N_{\mathrm{mode}}/2=128$ iterations, with the whole optimization running over
$3N_{\mathrm{mode}}=768$ iterations. Using a 6-point estimation in the
beginning, this strategy reduced the measurement time by factor 2.25 compared
to keeping the initial scheme till the end. Note that this technique is not
limited to non-linear feedback and can also be applied in the linear case. For
an increased precision, also higher order terms of the expansion of
$F(\delta)$ can be taken into account.