11institutetext: Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333
CA Leiden, the Netherlands
11email<EMAIL_ADDRESS>
# Statistical strong lensing. IV. Inferences with no individual source
redshifts

Alessandro Sonnenfeld 11

###### Abstract

Context.  Strong lensing mass measurements require the knowledge of the
redshift of both the lens and the source galaxy. Traditionally, spectroscopic
redshifts are used for this purpose. Upcoming surveys, however, will lead to
the discovery of $\sim 10^{5}$ strong lenses, and it will be very difficult to
obtain spectroscopic redshifts for most of them. Photometric redshift
measurements will also be very challenging due to the blending between lens
and source light.

Aims.  The goal of this work is to demonstrate how to carry out an inference
of the structural properties of the galaxy population from the analysis of a
set of strong lenses with no individual source redshift measurements, and to
assess the loss in precision compared to the case in which spectroscopic
redshifts are available.

Methods.  Building on the formalism introduced in Paper III, I developed a
method that allows a statistical strong lensing inference to be carried out
while marginalising over the source redshifts. This method, which relies on
the knowledge of the properties of the unlensed background source population
and of the selection function of the survey, generalises an approach known as
photogeometric redshift, originally introduced by the Strong Lensing Legacy
Survey collaboration. I tested the method on simulated data consisting of a
subset of 137 strong lenses that is complete above a cut in observational
space.

Results.  The method recovers the properties of the galaxy population with a
precision that is comparable to that attainable in the case in which
individual source redshifts are known.

Conclusions.  The photogeometric redshift method is a viable approach for the
analysis of large sets of strong lenses provided that the background source
population properties and lens selection function are well known.

###### Key Words.:

Gravitational lensing: strong – Galaxies: fundamental parameters

##  1 Introduction

Strong gravitational lensing allows for some of the most precise mass
measurements in astrophysics. The angular separation between strongly lensed
images depends on the mass distribution of the lens, as well as on the
distances between the observer, the lens, and the background source. When the
redshifts of the lens and source galaxies are known, by modelling the angular
structure of the strongly lensed images it is possible to obtain measurements
of the total projected mass enclosed within the images with a precision of a
few percent, as that is the typical uncertainty associated with the lens
modelling step (Bolton et al., 2008).

Traditionally, lens and source redshifts are obtained via spectroscopy.
Spectroscopic follow-up observations, however, can be expensive in terms of
exposure time. In the case of the Strong Lensing Legacy Survey (SL2S; Ruff et
al., 2011; Gavazzi et al., 2012), for example, 1-hour-long observations on
$8$-metre-class telescopes were needed to obtain the redshift of the
background sources. Moreover, most strongly lensed sources span the redshift
range $1-3$ (Sonnenfeld et al., 2013). A substantial fraction of them do not
have any strong emission lines in the optical region of the spectrum, and as a
result they go undetected when observed with optical spectroscopy (Ruff et
al., 2011; Sonnenfeld et al., 2013). Near-infrared observations usually
guarantee higher chances of success (Sonnenfeld et al., 2015, 2019), but very
few of the existing spectrographs have both the broad wavelength coverage
needed to maximise the chances of detecting an emission line and the high
spectral resolution required to minimise the contamination from atmospheric
emission.

Upcoming surveys such as Euclid111https://www.euclid-ec.org/, the Legacy
Survey of Space and Time (LSST222https://www.lsst.org/) and the Chinese Space
Station Telescope (CSST) will lead to the discovery of $\sim 10^{5}$ new
strong lenses (Collett, 2015). These lensing data have the potential to enable
unique measurements of the mass structure of galaxies. Sonnenfeld & Cautun
(2021 , hereafter Paper I) showed how with a sample of $10^{3}$ strong lenses
it is possible to calibrate stellar mass-to-light ratio measurements of
massive galaxies with $0.03$ dex precision. The analysis of Paper I, however,
was based on the assumptions that the redshifts of all lens galaxies and
background sources were known from spectroscopy. It is unclear whether such
spectroscopic data will become available in the near future.

The standard alternative to spectroscopy is photometric redshift (photo-z)
measurements. The accuracy of photo-zs depends critically on that of the
photometric measurements on which they are based. Most lenses are massive
elliptical galaxies, which are bright and with an easily characterised
spectral energy distribution. This makes it possible to obtain very precise
photometric redshift measurements for them (see for example Vakili et al.,
2019). Lensed sources, however, are fainter. Their images have typical
$i-$band magnitudes around $24-25$ and are usually found $1-2^{\prime\prime}$
away from the lens centre (Sonnenfeld et al., 2019). Consequently, they are
often blended with the light from the lens, which is usually brighter at all
but the bluest wavelengths. This makes it difficult to determine the colours
of lensed sources. As a result, photometric redshifts have so far found very
little use in galaxy-scale strong lensing, with the exception of cases in
which high-resolution imaging data in many bands are available (Sonnenfeld et
al., 2012).

When it is not possible to obtain individual redshift measurements, the only
remaining option is to rely on statistics: if the redshift distribution of the
population of background galaxies is known, then it is possible to guess the
source redshift of a lens, in a probabilistic sense. Since lensed sources
typically span a broad redshift range, this approach can no longer provide
percent mass measurements on individual systems. Nevertheless, high precision
can be gained with the statistical combination of a large number of lenses.

A similar approach is used in weak lensing analyses, where, for the
determination of the lensing critical surface mass density, the population
distribution of the source redshifts, ${\rm P}(z_{\mathrm{s}})$, is used
instead of individual photo-z measurements (see for example Hildebrandt et
al., 2017). In the case of strong lensing, however, there is an important
complication due to the selection function of the strong lens population. The
probability of a galaxy being strongly lensed by a foreground object varies
with its redshift, and so does the probability of it being detected by a
strong lens survey. As a result, the redshift distribution of strongly lensed
sources is very different from that of unlensed galaxies: for example, the
abundance of sources at a redshift close to that of the lenses is highly
suppressed. Therefore, it is not possible to directly use an estimate of ${\rm
P}(z_{\mathrm{s}})$ relative to the general background galaxy population as a
prior on the source redshift of a strong lens, as that would lead to a biased
measurement.

Ruff et al. (2011) tackled this problem in the following way in their study of
SL2S lenses. They assumed that galaxies have an isothermal density profile.
Then, they estimated a prior probability distribution for the lens deflection
strength, ${\rm P}(\sigma_{\mathrm{SIE}})$, of the lenses in the SL2S sample
by taking the observed stellar velocity dispersion function of galaxies and
re-weighting it to take into account the strong lensing cross-section and the
luminosity selection of the parent galaxy sample. The prior probability ${\rm
P}(\sigma_{\mathrm{SIE}})$ was then used, together with the observed Einstein
radius of each lens, to obtain a posterior probability distribution of the
source redshift. This method of obtaining the redshift of a strongly lensed
source is known as photogeometric redshift.

There are two problems with the Ruff et al. (2011) method. The first is that
it rests on the assumption of an isothermal density profile for the foreground
galaxy population. The density profile of galaxies, while on average close to
isothermal, varies with properties such as stellar mass, half-light radius,
and redshift (Sonnenfeld et al., 2013). Therefore, the isothermal assumption
introduces a bias. While this bias is probably not significant for the Ruff et
al. (2011) analysis, it is set to become an important source of systematic
errors once precision improves. More generally, if the purpose of the lensing
observations is to determine the structural properties of the galaxy
population, then the correct approach is to infer these properties jointly
with the redshifts of the background sources, instead of keeping them fixed.
The second problem is that Ruff et al. (2011) did not account for the
selection associated with the lens detection phase, which, in general, varies
with source redshift. Both of these problems can be overcome with the
formalism introduced in Sonnenfeld (2022 , Paper III from here on).

The purpose of this paper is to show how to carry out a population study of a
sample of strong lenses with no individual source redshift measurements, and
to quantify the loss in precision compared to the case in which spectroscopic
redshifts are available. The method that I propose relies on a good knowledge
of the properties of the unlensed background galaxy population, as well as the
selection function of the strong lensing survey. As with Paper III, this is
more easily done when the sample is complete above a well-defined cut in
observational space. I tested the method on a simulated population of
foreground galaxies and background sources, from which a complete subset of
strong lenses was drawn.

The structure of this work is the following. In Sect. 2 I introduce the
theoretical formalism on which this work is based. In Sect. 3 I describe the
simulations used to test the analysis method. In Sect. 4 I describe the model
that I used to fit the simulated data. In Sect. 5 I show the results of the
experiment. I discuss the results and draw conclusions in Sect. 6.

Throughout the paper I assume a flat $\Lambda$ cold dark matter cosmology with
$\Omega_{M}=0.3$ and $H_{0}=70\,\rm{km}\,\rm{s}^{-1}\,\rm{Mpc}^{-1}$. Masses
and luminosities are in solar units. The Python code used for the simulation
and analysis of the lens sample can be found on
GitHub333https://github.com/astrosonnen/strong_lensing_tools.

##  2 Photogeometric redshifts

This section explains (1) how information from strong lensing data, together
with prior knowledge of the background source population, can be used to
constrain the redshift of a lensed source, and (2) how the properties of a
population of galaxies can be inferred by statistically combining data from a
subset of lenses, even when their source redshifts are not known. The
following derivation builds on the formalism developed in Paper III.

The population of strong lenses from a survey with selection criterion $S$ can
be described with the following probability distribution:

${\rm
P}_{\mathrm{{SL}}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|S)\propto{\rm
P}_{\mathrm{g}}(\boldsymbol{\psi}_{\mathrm{g}}){\rm
P}_{\mathrm{s}}(\boldsymbol{\psi}_{\mathrm{s}}){\rm
P}_{\mathrm{det}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|S),$
(1)

where $\boldsymbol{\psi}_{\mathrm{g}}$ and $\boldsymbol{\psi}_{\mathrm{s}}$
are parameters describing the foreground galaxy and background source, ${\rm
P}_{\mathrm{g}}$ and ${\rm P}_{\mathrm{s}}$ are the galaxy and source
population distribution before lensing, and ${\rm P}_{\mathrm{det}}$ is the
probability of a lens-source pair being detected as a strong lens by the
survey. The source redshift, $z_{\mathrm{s}}$, is one of the parameters
describing the source. The source parameter array,
$\boldsymbol{\psi}_{\mathrm{s}}$, can then be written as

$\boldsymbol{\psi}_{\mathrm{s}}=(\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})},z_{\mathrm{s}}),$
(2)

where $\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})}$ includes all of the
source parameters (for instance, the magnitude, the position, etc.) except the
redshift.

Given a strong lens with observed data $\mathbf{d}$, the marginal posterior
probability of the source redshift is

${\rm P}(z_{\mathrm{s}}|\mathbf{d})\propto\int
d\boldsymbol{\psi}_{\mathrm{g}}d\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})}{\rm
P}_{\mathrm{{SL}}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})},z_{\mathrm{s}}|S){\rm
P}(\mathbf{d}|\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})},z_{\mathrm{s}}).$
(3)

If the strong lens population distribution is known exactly from prior
information, then one can use that knowledge to constrain ${\rm
P}(z_{\mathrm{s}}|\mathbf{d})$ via the above equation. This is the basic
principle of photogeometric redshift.

As the term photogeometric suggests, the constraints on the source redshift
consist of a photometric part, which is given by the knowledge, informed by
the lens model, of the unlensed surface brightness distribution of the source,
and a geometric part, which is directly related to the lens configuration. To
clarify this concept, it is useful to consider the following example.

We take a strong lens consisting of a foreground galaxy, described by
parameters $\boldsymbol{\psi}_{\mathrm{g}}$, and a background point source.
The data, $\mathbf{d}$, include the positions and apparent magnitudes of the
multiple images of the source. These data constrain the source redshift via
the factor ${\rm
P}(\mathbf{d}|\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})},z_{\mathrm{s}})$
in Eq. 3, which is the geometric term: modifying the source redshift varies
the predicted image positions and magnitudes as a result of the dependence of
the lensing deflection angle on the observer-lens-source geometry.

The photometric constraint on $z_{\mathrm{s}}$ is instead given by the factor
${\rm
P}_{\mathrm{{SL}}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})},z_{\mathrm{s}}|S)$
in Eq. 3. This can be written as

${\rm
P}_{\mathrm{{SL}}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})},z_{\mathrm{s}}|S)={\rm
P}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})}|S){\rm
P}(z_{\mathrm{s}}|\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})},S),$
(4)

where I used the rule of conditional probability. The first factor in the
right-hand side of the above equation is the prior probability of the lens and
source parameters (except the source redshift) being
$\boldsymbol{\psi}_{\mathrm{g}}$ and
$\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})}$ given the fact that the
lens is included in the survey with selection criterion $S$. The second
factor, ${\rm
P}(z_{\mathrm{s}}|\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})},S),$
is the prior probability of the source redshift given the remaining model
parameters and the selection criterion. One of the parameters entering
$\boldsymbol{\psi}_{\mathrm{s}}^{(-z_{\mathrm{s}})}$ is the intrinsic
(unlensed) source apparent magnitude, which can be used to further constrain
the source redshift given the knowledge of the distribution of background
sources in redshift-magnitude space. For example, if the source is
intrinsically very bright, large values of the redshift will be disfavoured.

Source redshift information is rarely needed for its own sake: in most cases,
$z_{\mathrm{s}}$ is required to constrain lens models. In other words, the
factor ${\rm P}_{\mathrm{g}}$, which enters Eq. 3 via the expansion of Eq. 1,
is uncertain, and the source redshift is one of the ingredients needed to
constrain it. Ruff et al. (2011) solved the problem of determining the
structural parameters of lenses with no source redshifts as follows. They used
information on ${\rm P}_{\mathrm{g}}$ that was independent of their strong
lens survey, the SL2S, together with the assumption of an isothermal density
profile for the galaxies, to infer the $z_{\mathrm{s}}$ of sources with no
spectroscopic redshifts. Then, they combined these estimates of
$z_{\mathrm{s}}$ with the observed Einstein radii and stellar velocity
dispersion to infer the evolution of the density profile of the lenses,
effectively obtaining an updated estimate of ${\rm P}_{\mathrm{g}}$. Their
procedure ensured that no data were used more than once, but it is internally
inconsistent because their inferred density profiles of the lenses were
different from isothermal.

The correct approach, when the main goal is to constrain the galaxy population
properties, is to infer ${\rm P}_{\mathrm{g}}$ while marginalising over the
source redshifts. If we summarise the galaxy distribution with a set of
parameters, $\boldsymbol{\eta}$, then we can obtain the posterior probability
distribution of $\boldsymbol{\eta}$ given the data as

$\displaystyle{\rm P}(\boldsymbol{\eta}|\mathbf{d})\propto{\rm
P}(\boldsymbol{\eta})\prod_{i}^{N_{\mathrm{lens}}}\int
d\boldsymbol{\psi}_{\mathrm{g},i}d\boldsymbol{\psi}_{\mathrm{s},i}^{(-z_{\mathrm{s}})}$
$\displaystyle dz_{\mathrm{s},i}{\rm
P}(\mathbf{d}_{i}|\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s},i}^{(-z_{\mathrm{s}})},z_{\mathrm{s}})\times$
$\displaystyle\times{\rm
P}_{\mathrm{{SL}}}(\boldsymbol{\psi}_{\mathrm{g},i},\boldsymbol{\psi}_{\mathrm{s},i}^{(-z_{\mathrm{s}})},z_{\mathrm{s},i}|\boldsymbol{\eta},S),$
(5)

where the product is taken over a sample of $N_{\mathrm{lens}}$ strong lenses.
This is the approach adopted in this paper.

The integrands in the right-hand side of Eq. 2 contain the expression for the
photogeometric redshift given by Eq. 3. In this case, however, the probability
distribution of the lens parameters, ${\rm P}_{\mathrm{{SL}}}$, is not fixed
but is a function of the galaxy population parameters, $\boldsymbol{\eta}$. If
a new lens is observed from a survey with the same characteristics, then its
photogeometric redshift can be estimated by marginalising over all possible
values of $\boldsymbol{\eta}$ given our knowledge gathered from the analysis
of the first $N_{\mathrm{lens}}$ lenses:

$\displaystyle{\rm
P}(z_{\mathrm{s},N_{\mathrm{lens}}+1}|\mathbf{d}_{N_{\mathrm{lens}}+1})\propto\int
d\boldsymbol{\eta}{\rm P}(\boldsymbol{\eta}|\mathbf{d})\int
d\boldsymbol{\psi}_{\mathrm{g},N_{\mathrm{lens}}+1}d\boldsymbol{\psi}_{\mathrm{s},N_{\mathrm{lens}}+1}^{(-z_{\mathrm{s}})}$
$\displaystyle{\rm
P}(\mathbf{d}_{N_{\mathrm{lens}}+1}|\boldsymbol{\psi}_{\mathrm{g},N_{\mathrm{lens}}+1},\boldsymbol{\psi}_{\mathrm{s},N_{\mathrm{lens}}+1}^{(-z_{\mathrm{s}})}){\rm
P}_{\mathrm{{SL}}}(\boldsymbol{\psi}_{\mathrm{g},N_{\mathrm{lens}}+1},\boldsymbol{\psi}_{\mathrm{s},N_{\mathrm{lens}}+1}^{(-z_{\mathrm{s}})}|\boldsymbol{\eta},S).$
(6)

##  3 Simulations

In order to demonstrate how the formalism introduced in the previous section
can be applied to a practical case, I generated a simulated sample of strong
lenses with the following general characteristics. As in the previous papers
of this series, I assumed foreground galaxies to have an axisymmetric mass
distribution and background sources to be point-like. Each galaxy consists of
the sum of two concentric components: a dark matter halo and a stellar bulge.
I defined the galaxy population among which strong lenses are searched for by
applying a luminosity cut. The galaxy sample is assumed to be complete above
this cut. Lens-source systems that produce at least two images with observed
magnitude brighter than a given threshold are listed as strong lenses detected
by the survey. The survey is assumed to be complete above this limit. In the
next few sections, I provide the details of how this sample was created.

###  3.1 Foreground galaxy population

I generated $10^{4}$ galaxies, with luminosity, $L$, drawn from the following
probability distribution

${\rm
P}(L)\sim\left(\frac{L}{L^{*}}\right)^{\alpha_{L}}\exp{\left\\{-\frac{L}{L^{*}}\right\\}},$
(7)

with $\log{L^{*}}=10.8$, $\alpha_{L}=-1.0$, and $\log{L}>10.5$. Given the
luminosity of each galaxy, I obtained its stellar mass by drawing the base-10
logarithm of the stellar mass-to-light ratio, $\Upsilon_{*}$, from the
following distribution:

$\log{\Upsilon_{*}}\sim\mathcal{N}(0.5,0.1^{2}),$ (8)

where the notation $\mathcal{N}(\mu,\sigma^{2})$ indicates a Gaussian
distribution with mean $\mu$ and variance $\sigma^{2}$. The combination of the
luminosity function, Eq. 7, and the mass-to-light ratio distribution of Eq. 8
produces a stellar mass function that is approximately that of quiescent
galaxies at $0.2<z<0.5$, as measured by Muzzin et al. (2013) under the
assumption of a Chabrier (2003) stellar initial mass function but shifted
towards higher masses by $\sim 0.1$ dex.

I described the mass distribution of the stellar component as an axisymmetric
de Vaucouleurs profile, with half-light radius drawn from the following
stellar mass-size relation:

$\log{R_{e}}\sim\mathcal{N}(1.0+0.8(\log{M_{*}}-11.5),0.15^{2}).$ (9)

Subsequently, I assigned a dark matter halo mass, $M_{h}$, to each galaxy. I
defined $M_{h}$ as the mass enclosed within a sphere with average density
equal to $200$ times the critical density of the Universe. I drew the value of
$\log{M_{h}}$ from the following Gaussian distribution:

$\log{M_{h}}\sim\mathcal{N}(13.0+1.5(\log{M_{*}}-11.5),0.2^{2}).$ (10)

The mean of this Gaussian scales with a power $1.5$ of the stellar mass,
$M_{*}$. I assumed a spherical Navarro, Frenk, and White (NFW; Navarro et al.,
1997) profile for the dark matter distribution, with a fixed ratio between the
virial radius and the scale radius, corresponding to a concentration of $5$:

$\frac{r_{200}}{r_{s}}=5,$ (11)

where $r_{200}$ is the radius of the sphere enclosing a dark matter mass equal
to $M_{h}$. For simplicity, the foreground galaxies are assumed to be all at
the same redshift, $z_{\mathrm{g}}=0.4$. The distribution in $\log{L}$,
$\log{M_{h}}$, $\log{\Upsilon}$, and $\log{R_{e}}$ of the sample of foreground
galaxies is shown as solid black contours in Fig. 1.

###  3.2 Background source population

I described each background source by means of its $i-$band apparent
magnitude, $m_{\mathrm{s}}$, and its redshift, $z_{\mathrm{s}}$. I first
assumed that their UV absolute magnitude, $M$, is drawn from the following
luminosity function:

$N(M)=0.4\ln{10}\phi^{*}\left(10^{-0.4(M-M^{*})}\right)^{\alpha_{s}+1}\exp{\left\\{-10^{-0.4(M-M^{*})}\right\\}},$
(12)

with $\phi^{*}=7.02\times 10^{-3}$ Mpc-3 mag-1, $M^{*}=-19.68$, and
$\alpha_{s}=-1.32$. These values correspond to the rest-frame $1500\AA$
luminosity function of $z\approx 2$ galaxies, as measured by Parsa et al.
(2016).

Then, I assumed that all of the sources have the same spectral energy
distribution, given by the spectral template ‘starb1’ from the Kinney-Calzetti
Spectral Atlas of Galaxies (Calzetti et al., 1994; Kinney et al., 1996). For
ease of computation, I only considered galaxies with $m_{\mathrm{s}}<27$ and
$z_{\mathrm{g}}+0.1<z_{\mathrm{s}}<4$. The resulting distribution in
$z_{\mathrm{s}}$ and $m_{\mathrm{s}}$ is shown as solid black lines in Fig. 1.
Once integrated over $z_{\mathrm{s}}$, the background source population has an
average number density of $n_{\mathrm{bkg}}=162\,\mathrm{arcmin}^{-2}$.

###  3.3 Strong lens definition

Given each foreground galaxy, I randomly placed background sources behind it,
assuming a uniform spatial distribution with average projected number density
$n_{\mathrm{bkg}}$. Then, I solved the lens equation to determine their
observed image positions. For simplicity, I assumed the background sources to
be massless so that the only source of lensing deflection was that of the
foreground galaxy. Axisymmetric lenses of the kind considered in this
experiment can produce up to three images of the same source. I focused on the
two brightest images, at positions $\theta_{1}$ and $\theta_{2}$. With the
label $1$ I indicate the main image, that is, the one outside of the
tangential critical curve, which is usually the brightest. Image $2$, instead,
is located between the radial and the tangential critical curve, opposite to
image $1$ with respect to the lens centre (which is the same notation used in
all of the papers of this series).

Given the intrinsic apparent magnitude, $m_{\mathrm{s}}$, of a strongly lensed
source and the lensing magnification at the two image positions, I computed
the magnitudes of the two lensed images, $m_{1}$ and $m_{2}$. Then, I added a
Gaussian observational noise with dispersion $\epsilon_{m}=0.1$ mag to obtain
their corresponding observed values:

$\displaystyle m_{1}^{\mathrm{obs}}$
$\displaystyle\sim\mathcal{N}(m_{1},\epsilon_{m}^{2})$ (13) $\displaystyle
m_{2}^{\mathrm{obs}}$ $\displaystyle\sim\mathcal{N}(m_{2},\epsilon_{m}^{2}).$
(14)

I defined as a strong lens any galaxy with a strongly lensed source for which
the observed magnitude of image $2$ is brighter than $26$. The lens selection
criterion is then

$S=\left\\{m_{2}^{\mathrm{obs}}<26\right\\}.$ (15)

This condition mimics a selection based on the detection of a counter-image to
the main arc, which usually allows the lensing nature of a lens candidate to
be determined with high confidence. I assumed that all such lenses are
included in the survey, which is therefore complete above the observational
cut given by Eq. 15. All of the other galaxies are treated as non-detections.

With this definition, the simulation produced $137$ strong lenses. Two of the
lenses have more than one strongly lensed source: in those cases, I only
considered one lensed source each, chosen at random. The distribution in
foreground galaxy properties, background source properties, and Einstein
radius of the strong lenses is shown in Fig. 1. Compared to the general
population, strong lens galaxies tend to be more luminous and, consequently,
to have higher halo masses. Additionally, they are more compact at fixed
luminosity (that is, their half-light radius is on average smaller) and have a
higher stellar mass-to-light ratio. These last two trends are quantified in
Table 1.

The distribution of strongly lensed sources differs from the general
population of background galaxies in two ways. First, the low-redshift part of
the distribution is suppressed. This is because the strong lensing cross-
section decreases with source redshift at fixed lens redshift. Second, the
abundance of faint galaxies ($m_{\mathrm{s}}>26$) is greatly reduced. This is
a consequence of the detection limit of Eq. 15.

Figure 1:  Simulated sample. Solid lines show the distribution in the
luminosity, halo mass, stellar mass-to-light ratio, and half-light radius of
the general population of foreground galaxies as well as the distribution in
redshift and $i-$band magnitude of the general population of background
sources. Contour levels correspond to 68% and 95% enclosed probability. Pink
points show the values of the above quantities, as well as of the Einstein
radius, of the strong lenses. Histograms are normalised to one.

##  4 Model

I fitted the simulated strong lensing data with a model that describes the
joint distribution of the galaxy population, the background source population,
and the strong lens detection efficiency of the survey, using the same method
introduced in Paper III. Each aspect of the model is described in the
following sections.

###  4.1 Foreground galaxies

I described each galaxy with the sum of a stellar component, consisting of a
de Vaucouleurs profile with known luminosity and half-light radius and unknown
mass-to-light ratio, and a dark matter halo, consisting of an NFW profile.
Moreover, I assumed that the concentration parameter of the NFW profile, Eq.
11, is known exactly from prior information. This is not a realistic
assumption, but it is a useful one: it simplifies calculations while still
allowing the main goal of this paper to be pursued, which is to demonstrate
the accuracy of the inference method itself.

Each galaxy is therefore described by four parameters: the luminosity, $L$,
the stellar mass-to-light ratio, $\Upsilon_{*}$, the half-light radius,
$R_{e}$, and the halo mass, $M_{h}$. The galaxy parameter set is therefore

$\boldsymbol{\psi}_{\mathrm{g}}\equiv\\{L,\Upsilon_{*},R_{e},M_{h}\\}.$ (16)

I assumed the following probability distribution for
$\boldsymbol{\psi}_{\mathrm{g}}$:

${\rm
P}_{\mathrm{g}}(\boldsymbol{\psi}_{\mathrm{g}}|\boldsymbol{\eta})=\mathcal{L}(L)\mathcal{U}(\log{\Upsilon_{*}})\mathcal{R}(\log{R_{e}})\mathcal{H}(\log{M_{h}}).$
(17)

The factor $\mathcal{L}$ is the luminosity distribution of the sample, which I
fixed to the true distribution used to generate the mock. The factor
$\mathcal{U}$ is the distribution in the stellar mass-to-light ratio, which I
assumed to be a Gaussian in $\log{\Upsilon_{*}}$,

$\mathcal{U}(\log{\Upsilon_{*}})\sim\mathcal{N}(\mu_{\Upsilon},\sigma_{\Upsilon}^{2}),$
(18)

with mean $\mu_{\Upsilon}$ and dispersion $\sigma_{\Upsilon}$.

The factor $\mathcal{R}$ is the distribution in half-light radius, which I
modelled as the following Gaussian in $\log{R_{e}}$:

$\mathcal{R}(\log{R_{e}})\sim\mathcal{R}(\mu_{\mathrm{R},0}+\beta_{\mathrm{R}}(\log{M_{*}}-\log{M_{*}^{p}}),\sigma_{\mathrm{R}}^{2}),$
(19)

with

$\log{M_{*}^{p}}=11.0+\mu_{\Upsilon}.$ (20)

In words, the mean $\log{R_{e}}$ scales with a power $\beta_{R}$ of the
stellar mass, and the pivot point of the stellar mass-size relation,
$M_{*}^{p}$, is the average stellar mass of a galaxy with luminosity
$\log{L}=11$. The parameter $\mu_{\mathrm{R},0}$, then, is the average
$\log{R_{e}}$ of a galaxy with luminosity $\log{L}=11$. This quantity can be
measured with high precision on the sample of galaxies from which lenses are
drawn; therefore, I assume that it is known exactly.

Finally, the factor $\mathcal{H}$ is the distribution in halo mass. Similarly
to the $\mathcal{R}$ factor, this is modelled as

$\mathcal{H}(\log{M_{h}})\sim\mathcal{N}(\mu_{h,0}+\beta_{h}(\log{M_{*}}-\log{M_{*}^{p}}),\sigma_{h}^{2}).$
(21)

The distribution ${\rm
P}_{\mathrm{g}}(\boldsymbol{\psi}_{\mathrm{g}}|\boldsymbol{\eta})$ is a
function of the following set of free parameters:

$\boldsymbol{\eta}\equiv\left\\{\mu_{\Upsilon},\sigma_{\Upsilon},\beta_{\mathrm{R}},\sigma_{\mathrm{R}},\mu_{\mathrm{h},0},\beta_{\mathrm{h}},\sigma_{\mathrm{h}}\right\\}.$
(22)

The goal is to infer $\boldsymbol{\eta}$ given the data. The distribution
${\rm P}_{\mathrm{g}}(\boldsymbol{\psi}_{\mathrm{g}}|\boldsymbol{\eta})$
reduces to the true distribution used to generate the mock for

$\boldsymbol{\eta}^{(\mathrm{true})}=\left\\{0.5,0.1,0.8,0.15,13.0,1.5,0.2\right\\}.$
(23)

###  4.2 Background sources

Sources are point-like, and therefore they can be fully described with their
apparent magnitude, $m_{\mathrm{s}}$, their redshift, $z_{\mathrm{s}}$, and
their position relative to the lens centre, $\beta$:

$\boldsymbol{\psi}_{\mathrm{s}}\equiv\left\\{m_{\mathrm{s}},z_{\mathrm{s}},\beta\right\\}.$
(24)

I assumed that the distribution in $\boldsymbol{\psi}_{\mathrm{s}}$ of the
background galaxy population is known exactly from prior observations. I
discuss possible strategies to drop this assumption in Sect. 6.

The probability distribution ${\rm
P}_{\mathrm{s}}(\boldsymbol{\psi}_{\mathrm{s}})$ is therefore fixed to the one
used to generate the mock (shown in black in Fig. 1). This can be written as

${\rm
P}_{\mathrm{s}}(\boldsymbol{\psi}_{\mathrm{s}})=U_{\mathbb{R}^{2}}(\boldsymbol{\beta}){\rm
P}(z_{\mathrm{s}}){\rm P}(m_{\mathrm{s}}|z_{\mathrm{s}}).$ (25)

In the above equation, $U_{\mathbb{R}^{2}}$ is a uniform distribution on the
sky, ${\rm P}(m_{\mathrm{s}}|z_{\mathrm{s}})$ is

${\rm
P}(m_{\mathrm{s}}|z_{\mathrm{s}})\propto\left(10^{-0.4\left(m_{\mathrm{s}}-m_{\mathrm{s}}^{*}(z_{\mathrm{s}})\right)}\right)^{\alpha_{s}+1}\exp{\left\\{-10^{-0.4\left(m_{\mathrm{s}}-m_{\mathrm{s}}^{*}(z_{\mathrm{s}})\right)}\right\\}},$
(26)

where $m_{\mathrm{s}}^{*}(z_{\mathrm{s}})$ is obtained from $M^{*}$ via the
same spectral template used to generate the mock, and ${\rm
P}(z_{\mathrm{s}})$ is

${\rm P}(z_{\mathrm{s}})\propto\frac{dV_{C}}{dz}\frac{dn_{\mathrm{bkg}}}{dz},$
(27)

which is the product of the derivatives with respect to the redshift of the
comoving volume and the surface density of background sources.

###  4.3 Lens detection efficiency

As Sect. 3.3 explains, the lens sample is complete above the cut of Eq. 15.
This implies that the detection efficiency of a strong lens, expressed in
terms of model quantities, is given by

${\rm
P}_{\mathrm{det}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|S)=\frac{1}{2}\left[1-{\mathrm{erf}}{\left(\frac{m_{2}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}})-26}{\sqrt{2}\epsilon_{m}}\right)}\right].$
(28)

###  4.4 Inferring the galaxy population parameters

I fitted the model described in this section to the data relative to the $137$
strong lenses. For each lens, I assumed that the following data are available:
(1) the positions of the two images of the lensed source,
$\theta_{1}^{\mathrm{obs}}$ and $\theta_{2}^{\mathrm{obs}}$, measured exactly;
(2) the half-light radius of the galaxy, $R_{e}^{(\mathrm{obs})}$, measured
exactly; (3) the luminosity of the lens galaxy, $L^{(\mathrm{obs})}$, measured
exactly; and (4) the magnitudes of the two images, $m_{1}^{\mathrm{obs}}$ and
$m_{2}^{\mathrm{obs}}$, measured with uncertainty $\epsilon_{m}$.

In principle, the remaining galaxies in the sample (the non-detections) could
be used as additional constraints, following the method introduced in Paper
III. However, Paper III also showed that, when relative magnification
information is available (that is, when the magnitudes of the lensed images
are used as constraints), the benefit of adding the non-detections is
negligible.

Following Paper III, the posterior probability distribution of the galaxy
population parameters, given the data of the strong lenses, is given by

${\rm P}(\boldsymbol{\eta}|\mathbf{d})\propto{\rm
P}(\boldsymbol{\eta})\prod_{i}{\rm P}(\mathbf{d}_{i}|\boldsymbol{\eta}),$ (29)

where the product is taken over the strong lenses, and $\mathbf{d}_{i}$ is the
data relative to the $i-$th lens. Each factor ${\rm
P}(\mathbf{d}_{i}|\boldsymbol{\eta})$ can be written as the following
integral:

${\rm P}(\mathbf{d}_{i}|\boldsymbol{\eta})=\int
d\boldsymbol{\psi}_{\mathrm{g}}d\boldsymbol{\psi}_{\mathrm{s}}{\rm
P}(\mathbf{d}_{i}|\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}){\rm
P}_{\mathrm{{SL}}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|\boldsymbol{\eta},S),$
(30)

where

${\rm
P}_{\mathrm{{SL}}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|\boldsymbol{\eta},S)=A(\boldsymbol{\eta}){\rm
P}_{\mathrm{g}}(\boldsymbol{\psi}_{\mathrm{g}}|\boldsymbol{\eta}){\rm
P}_{\mathrm{s}}(\boldsymbol{\psi}_{\mathrm{s}}){\rm
P}_{\mathrm{det}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|S),$
(31)

and $A(\boldsymbol{\eta})$ is a constant that ensures that ${\rm
P}_{\mathrm{{SL}}}$ is normalised to unity. The likelihood of observing the
data given the lens and source parameters, ${\rm
P}(\mathbf{d}_{i}|\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}})$,
is

$\displaystyle{\rm
P}(\mathbf{d}_{i}|\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}})=$
$\displaystyle\delta\left(R_{e}-R_{e}^{(\mathrm{obs})}\right)\delta\left(L-L^{(\mathrm{obs})}\right)\delta\left(\theta_{1}(\boldsymbol{\psi}_{\mathrm{g}},\beta)-\theta_{1}^{\mathrm{obs}}\right)$
$\displaystyle\delta\left(\theta_{2}(\boldsymbol{\psi}_{\mathrm{g}},\beta)-\theta_{2}^{\mathrm{obs}}\right){\rm
P}(m_{1}^{\mathrm{obs}},m_{2}^{\mathrm{obs}}|\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}),$
(32)

where ${\rm
P}(m_{1}^{\mathrm{obs}},m_{2}^{\mathrm{obs}}|\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}})$
is the product of two Gaussians, with means $m_{1}^{\mathrm{obs}}$ and
$m_{2}^{\mathrm{obs}}$ and dispersion $\epsilon_{m}$.

Equation 30 is a seven-dimensional integral. The integrals over the galaxy
half-light radius and luminosity are trivial, in virtue of the Dirac delta
functions of Eq. 4.4. Similarly, the fact that the two image positions are
known exactly allows the integration to be performed over two lens parameters
analytically. After integrating over the source position $\beta$ and
$\log{\Upsilon_{*}}$, Eq. 4.4 becomes

$\displaystyle{\rm P}(\mathbf{d}_{i}|\boldsymbol{\eta})=\int$ $\displaystyle
dz_{\mathrm{s}}dm_{\mathrm{s}}d\log{M_{h}}{\rm
P}(m_{1}^{\mathrm{obs}},m_{2}^{\mathrm{obs}}|M_{h},L\Upsilon_{*}^{*},R_{e},\beta^{*})$
$\displaystyle\times A(\boldsymbol{\eta}){\rm
P}(L,R_{e},M_{h},\Upsilon_{*}^{*}|\boldsymbol{\eta}){\rm
P}_{\mathrm{s}}(z_{\mathrm{s}},m_{\mathrm{s}},\beta^{*})$
$\displaystyle\times{\rm
P}_{\mathrm{det}}(\boldsymbol{\psi}_{\mathrm{g}},z_{\mathrm{s}},m_{\mathrm{s}},\beta^{*}|S)\left\lvert\mathrm{det}J\right\rvert_{\left(\Upsilon_{*},\beta\right)=\left(\Upsilon_{*}^{*},\beta^{*}\right)},$
(33)

where I dropped the subscript $i$ and set $L=L^{(\mathrm{obs})}$ and
$R_{e}=R_{e}^{(\mathrm{obs})}$ for ease of notation. In the above equation,
$\Upsilon_{*}^{*}$ and $\beta^{*}$ are the values of the stellar mass-to-light
ratio and source position needed to reproduce the two image positions, while
the last factor in the integrand is the Jacobian determinant of the following
variable change:

$(\log{\Upsilon_{*}},\beta)\rightarrow(\theta_{1},\theta_{2}).$ (34)

I sampled the posterior probability distribution of $\boldsymbol{\eta}$ via a
Markov chain Monte Carlo. I assumed a uniform prior on all of the parameters,
over the range listed in Table 1. At each draw of $\boldsymbol{\eta}$, I
computed both the normalisation constant, $A(\boldsymbol{\eta})$, and the
integrals of Eq. 4.4 via Monte Carlo integration with importance sampling.

###  4.5 The redshift-lens mass degeneracy

Before showing the results of the inference over the entire lens population,
it is useful to consider the problem of the degeneracy between the source
redshift and the lens mass parameters in the case of a single lens. The lens
model used for the fit has five degrees of freedom: the stellar and halo mass
of the lens galaxy and the magnitude, redshift, and position of the background
source. If the magnitudes of the two images and the redshift of the source
were known exactly (in addition to the image positions), then the model would
be fully constrained, and the result of the fit would be equal to the truth.
In the case of unknown source redshift, however, the uncertainty on
$z_{\mathrm{s}}$ propagates to the other lens model parameters.

In order to understand the amplitude of this effect, I fitted a lens from the
simulated sample, assuming perfect measurements, while assuming a range of
values of the source redshift. Figure 2 shows the inferred values of
$\Upsilon_{*}$, $M_{h}$ as a function of $z_{\mathrm{s}}$.

Figure 2:  Stellar mass-to-light ratio (upper panel) and halo mass (lower
panel) obtained by fitting the lens model introduced in Sect. 3.1 to a mock
strong lens from the simulated sample as a function of the assumed source
redshift. The fit was done assuming perfect measurements of the two image
positions and magnitudes.

Both the inferred stellar mass-to-light ratio and the halo mass increase with
decreasing source redshift: this is because, as the source gets closer to the
lens, the value of the surface critical density for strong lensing becomes
larger, and therefore a higher mass is required to produce the deflection
angles needed to fit the image positions. More precisely, when varying the
source redshift from $z_{\mathrm{s}}=1$ to $z_{\mathrm{s}}=3$, the inferred
stellar mass-to-light ratio changes from $\log{\Upsilon_{*}}=0.55$ to
$\log{\Upsilon_{*}}=0.42$. While this is a relatively large uncertainty, the
combination of a large set of lenses allows this source of noise to be be
beaten down statistically, as the next section shows.

##  5 Results

###  5.1 Galaxy population inference

Figure 3 shows the posterior probability distribution of the parameters that
describe the halo mass and stellar mass-to-light ratio distribution of the
galaxy population. The median and 68% credible region of the marginal
posterior of all of the model parameters (including those describing the size
distribution) is reported in Table 1. The true values of all of the parameters
are accurately recovered. The inference on the average $\log{\Upsilon_{*}}$
parameter, in particular, is $\mu_{\Upsilon}=0.498\pm 0.018$. Moreover, as in
Paper III, the model is able to correct for the lensing selection bias: it
favours values of the parameters closer to those of the general galaxy
population and disfavours values measured on the lenses alone (these values
are shown in Fig. 3 and reported in Table 1), though with low significance.

Figure 3:  Posterior probability distribution of the parameters describing the
distribution in $M_{h}$ and $\Upsilon_{*}$ of the galaxy population, obtained
by fitting the simulated sample of 137 strong lenses with the model of Sect.
4. Solid black contours show the inference obtained with no individual source
redshift information. Red filled contours show the inference obtained with
known individual source redshifts. Dashed black lines mark the true values of
the parameters. Dotted green lines indicate the values of the parameters
obtained by fitting the model to the subsample of strong lenses.  Table 1:
Inference on the galaxy population model parameters.

Parameter | Truth | Lenses | Prior | No individual | Known individual | Description
---|---|---|---|---|---|---
|  |  |  | redshifts | redshifts |
$\mu_{\mathrm{h},0}$ | $13.00$ | $13.08$ | $U(12.5,13.5)$ | $12.96\pm 0.06$ | $12.96\pm 0.07$ | Mean $\log{M_{h}}$ at luminosity $\log{L}=11$ and
|  |  |  |  |  | average $\Upsilon_{*}$
$\beta_{\mathrm{h}}$ | $1.50$ | $1.54$ | $U(0.0,3.0)$ | $1.48\pm 0.17$ | $1.57\pm 0.18$ | Scaling of the mean $\log{M_{h}}$ with $\log{M_{*}}$
$\sigma_{\mathrm{h}}$ | $0.20$ | $0.20$ | $U(0.0,0.5)$ | $0.21\pm 0.05$ | $0.22\pm 0.05$ | Scatter in $\log{M_{h}}$ around the mean
$\mu_{\Upsilon}$ | $0.50$ | $0.52$ | $U(0.3,0.7)$ | $0.494\pm 0.019$ | $0.514\pm 0.015$ | Mean $\log{\Upsilon_{*}}$
$\sigma_{\Upsilon}$ | $0.10$ | $0.10$ | $U(0.0,0.3)$ | $0.115\pm 0.012$ | $0.099\pm 0.007$ | Scatter in $\log{\Upsilon_{*}}$ around the mean
$\beta_{\mathrm{R}}$ | $0.80$ | $0.70$ | $U(0.0,2.0)$ | $0.76\pm 0.05$ | $0.79\pm 0.05$ | Scaling of $\log{R_{e}}$ with $\log{M_{*}}$
$\sigma_{\mathrm{R}}$ | $0.15$ | $0.14$ | $U(0.0,0.5)$ | $0.164\pm 0.011$ | $0.161\pm 0.010$ | Scatter in $\log{R_{e}}$ around the mean

444 Column (2): True values of the population parameters. Column (3): Values
of the lens population parameters, obtained by fitting the model of Eq. 17 to
the sample of detected strong lenses. Column (4): Priors on the parameters.
Column (5): Median and 68% credible interval of the marginal posterior
probability distribution of each parameter. Column (6): Inference with
individual source redshift information.

In order to compare this result with a case in which source redshifts are
available, I repeated the analysis by fixing the $z_{\mathrm{s}}$ of each lens
to the true value. The resulting posterior probability distribution is also
shown in Fig. 3 and reported in Table 1. The uncertainties obtained when
$z_{\mathrm{s}}$ is known are comparable to the case with no individual
$z_{\mathrm{s}}$ information. Only the inference on the parameters that
describe the stellar mass-to-light ratio is slightly more precise. This means
that the primary factor that determines the uncertainty in this experiment is
the limited amount of information that the lensing data provide, and not the
lack of source redshift information. The uncertainties associated with the
individual source redshifts are effectively reduced by the statistical
combination of the lenses in the sample.

This result depends on the relative importance of the source redshift
uncertainty in the error budget. In turn, this can vary with the quality of
the lensing data, the flexibility of the model used to fit the data, the
intrinsic scatter in the lens and source population, and the sample size.
Exploring the sensitivity of the results to these various aspects is beyond
the purpose of this paper. We can consider, however, the case in which the
two-component galaxy model of Sect. 3.1 is replaced with a singular isothermal
sphere (SIS). An SIS is fully described by a single parameter that, in the
case in which the source redshift is known, can be constrained exactly with
the two image positions of the lensed source. In other words, uncertainties
associated with the individual lens models vanish when assuming an SIS. As a
result, the uncertainty associated with the source redshifts acquires a bigger
weight in the error budget, and therefore the difference in precision between
the known and unknown individual redshift cases becomes larger. Nevertheless,
massive galaxies are known to have a more complex mass structure than a simple
SIS (see e.g. Sonnenfeld et al., 2013): in all practical cases, when analysing
real samples of lenses, models with at least two degrees of freedom are to be
preferred.

###  5.2 Photogeometric redshifts on new lenses

Once the properties of the population of galaxies have been determined, it is
possible to use this knowledge to predict the source redshifts of individual
strong lenses. We can consider a new lens, discovered by a survey with
identical properties to the one on which the mock sample was based. The
marginal posterior probability of the source redshift of this lens is given by
Eq. 2. Figure 4 shows this posterior probability for an example lens, fitted
with the same model of Sect. 4 and using the posterior probability of the
galaxy population parameters, ${\rm P}(\boldsymbol{\eta}|\mathbf{d})$,
obtained earlier.

In order to illustrate the way in which the knowledge of the galaxy population
informs this quantity, Fig. 4 also shows the marginal posterior probability on
$z_{\mathrm{s}}$ obtained (1) assuming flat priors on the lens parameters
$\log{M_{h}}$ and $\log{\Upsilon_{*}}$ and (2) neglecting the lens detection
efficiency term, ${\rm P}_{\mathrm{det}}$. The two distributions are very
different: the information on the galaxy population properties and on the
selection function correctly disfavours solutions with source redshifts very
close to the lens. On the contrary, when this information is not taken into
account, the posterior probability in the source redshift stays similar to the
prior distribution of the background source population, which is highest at
the lowest redshifts.

Figure 4:  Photogeometric redshift of a strongly lensed source, assuming a
lens model from the family of Sect. 3.1. The blue histogram was obtained by
fitting the lensing data assuming a flat prior on $\log{M_{h}}$ and
$\log{\Upsilon_{*}}$ and the source population distribution of Sect. 3.2, with
no correction for the strong lensing selection. The orange histogram was
obtained assuming the full model, including the selection function correction,
and using the values of the galaxy population parameters, $\boldsymbol{\eta}$,
inferred in the experiment as a prior. The green histogram was obtained
neglecting observational errors and fixing the values of $\boldsymbol{\eta}$
to the truth. The dashed vertical line marks the true value of the source
redshift.

The width of the ${\rm P}(z_{\mathrm{s}}|\mathbf{d}_{N_{\mathrm{lens}}+1})$
distribution is the result of (1) uncertainties associated with the lens
model, (2) uncertainties on the galaxy population distribution parameters,
$\boldsymbol{\eta}$, and (3) the width of the ${\rm
P}_{\mathrm{{SL}}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|\boldsymbol{\eta})$
distribution. In order to isolate these different sources of uncertainty, I
repeated the analysis assuming perfect lensing measurements (that is, assuming
no errors on the magnitudes of the two images) and fixing $\boldsymbol{\eta}$
to the true value. The resulting posterior probability distribution, also
shown in Fig. 4, is very similar to the one obtained without these additional
assumptions. This means that the main sources of uncertainty on the
photogeometric redshift are the width of the galaxy and source population
distributions.

##  6 Discussion and summary

The analysis described in the previous sections relied on the assumption of a
perfect knowledge of the distribution of the background source population,
${\rm P}_{\mathrm{s}}(\boldsymbol{\psi}_{\mathrm{s}})$. As discussed in Paper
III, it can be challenging to acquire this information in practice, especially
if the selection criteria of a survey involve more aspects than just the
magnitudes and positions of the lensed images. For example, if the detection
probability is also a function of surface brightness, then the model
$\boldsymbol{\psi}_{\mathrm{s}}$ must also include parameters such as galaxy
size.

In principle, ${\rm P}_{\mathrm{s}}(\boldsymbol{\psi}_{\mathrm{s}})$ can be
determined by studying the population of unlensed galaxies. However, because
of lensing magnification, the appearance of strongly lensed sources, and
therefore their detectability, can vary with changes in
$\boldsymbol{\psi}_{\mathrm{s}}$ that occur below the resolution limit of the
observations with which the unlensed galaxy population is studied. A possible
solution, then, is to determine ${\rm
P}_{\mathrm{s}}(\boldsymbol{\psi}_{\mathrm{s}})$ along with the galaxy
population, adding flexibility to the model. This will introduce an additional
source of uncertainty. In that case, the redshift dependence of ${\rm
P}_{\mathrm{s}}(\boldsymbol{\psi}_{\mathrm{s}})$ could be constrained by
obtaining spectroscopic redshifts for a small subset of the background sources
in the strong lens sample. I leave the exploration of the case in which the
source population distribution term is no longer fixed for future study.

This paper revisited the concept of photogeometric redshift, addressing the
theoretical shortcomings present in the original formulation of Ruff et al.
(2011). The experiment carried out demonstrates that it is possible to use
strong lenses with no individual source redshift information to make accurate
inferences of the properties of the galaxy population. When dealing with a
large sample of lenses ($N_{\mathrm{lens}}>100$), the statistical noise from
the uncertainty on the individual source redshifts is greatly reduced. In the
case in which the properties of the background (unlensed) source population
are known exactly, its contribution to the error budget on the galaxy
population parameters is only secondary.

The key elements that are necessary for a successful application of this
method are (1) the knowledge of the lens selection function and (2) the
knowledge of the distribution in surface brightness profile-redshift space of
the background galaxy population. The problem of determining the former can be
greatly simplified by identifying a subset of the strong lens population that
is complete above an observational cut, as shown in Paper III. Measuring the
latter is potentially made difficult by finite resolution effects, the
importance of which depends on the specific method used to discover lenses.
Nevertheless, the distribution of the source properties can be inferred along
with the properties of the foreground galaxy population by extending the
model. In summary, the method introduced in this paper provides a viable
solution to the problem of analysing large samples of lenses when
spectroscopic follow-up is not possible.

###### Acknowledgements.

I thank Phil Marshall and Raphael Gavazzi for useful suggestions.

## References

  * Bolton et al. (2008) Bolton, A. S., Burles, S., Koopmans, L. V. E., et al. 2008, ApJ, 682, 964
  * Calzetti et al. (1994) Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Collett (2015) Collett, T. E. 2015, ApJ, 811, 20
  * Gavazzi et al. (2012) Gavazzi, R., Treu, T., Marshall, P. J., Brault, F., & Ruff, A. 2012, ApJ, 761, 170
  * Hildebrandt et al. (2017) Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454
  * Kinney et al. (1996) Kinney, A. L., Calzetti, D., Bohlin, R. C., et al. 1996, ApJ, 467, 38
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Parsa et al. (2016) Parsa, S., Dunlop, J. S., McLure, R. J., & Mortlock, A. 2016, MNRAS, 456, 3194
  * Ruff et al. (2011) Ruff, A. J., Gavazzi, R., Marshall, P. J., et al. 2011, ApJ, 727, 96
  * Sonnenfeld (2022) Sonnenfeld, A. 2022, A&A, 659, A132, (Paper III)
  * Sonnenfeld & Cautun (2021) Sonnenfeld, A. & Cautun, M. 2021, A&A, 651, A18, (Paper I)
  * Sonnenfeld et al. (2019) Sonnenfeld, A., Jaelani, A. T., Chan, J., et al. 2019, A&A, 630, A71
  * Sonnenfeld et al. (2012) Sonnenfeld, A., Treu, T., Gavazzi, R., et al. 2012, ApJ, 752, 163
  * Sonnenfeld et al. (2013) Sonnenfeld, A., Treu, T., Gavazzi, R., et al. 2013, ApJ, 777, 98
  * Sonnenfeld et al. (2015) Sonnenfeld, A., Treu, T., Marshall, P. J., et al. 2015, ApJ, 800, 94
  * Vakili et al. (2019) Vakili, M., Bilicki, M., Hoekstra, H., et al. 2019, MNRAS, 487, 3715


# Iterative PDE-constrained optimization for seismic full-waveform inversion

Mikhail Malovichko  Moscow Institute of Physics and Technology  Akbar
Orazbayev  Moscow Institute of Physics and Technology  Nikolay Khokhlov
Moscow Institute of Physics and Technology

###### Abstract

This paper presents a novel numerical method for the Newton seismic full-
waveform inversion (FWI). The method is based on the full-space approach,
where the state, adjoint state, and control variables are optimized
simultaneously. Each Newton step is formulated as a PDE-constrained
optimization problem, which is cast in the form of the Karush-Kuhn-Tucker
(KKT) system of linear algebraic equitations. The KKT system is solved
inexactly with a preconditioned Krylov solver. We introduced two
preconditioners: the one based on the block-triangular factorization and its
variant with an inexact block solver. The method was benchmarked against the
standard truncated Newton FWI scheme on a part of the Marmousi velocity model.
The algorithm demonstrated a considerable runtime reduction compared to the
standard FWI. Moreover, the presented approach has a great potential for
further acceleration. The central result of this paper is that it establishes
the feasibility of Newton-type optimization of the KKT system in application
to the seismic FWI.

##  1 Introduction

The seismic full-waveform inversion (FWI) is a coefficient inverse problem
aiming to estimate the subsurface distribution of material properties based on
recorded seismic data. The FWI has emerged almost four decades ago [1, 2, 3],
and has evolved into a powerful method for reconstruction of the subsurface
properties, see reviews [4, 5].

The FWI is formulated as a Tikhonov-like minimization procedure,

$\underset{\gamma}{\text{minimize}}\quad\left\lVert
d-\mathcal{F}(\gamma)\right\rVert^{2}+\mathcal{R}(\gamma),$ (1)

where $\gamma$ is a function of material parameters (say, the squared P-wave
slowness), $d$ is a vector of observed data, $\mathcal{F}$ is a nonlinear
forward-problem operator, and $\mathcal{R}$ is the stabilizing term. We are
interested in the frequency-domain acoustic FWI and, thus, operator
$\mathcal{F}$ implies the solution of Helmholtz’s equation. Many variants of
basic formulation (1) exist: $\mathcal{R}$ can take various forms or even be
omitted, the residual can be weighted, and so on.

Optimization (1) is usually performed by either quasi-Newton or Newton
methods. The quasi-Newton methods dominate industrial applications due to
their relative numerical efficiency [6, 7, 8, 9, 10] (among many others).
These methods account for nonlinearity partially, usually by approximating the
Hessian with a diagonal matrix. Typical representatives of this class are the
nonlinear conjugate gradient (NLCG) method, where a diagonal matrix scales the
gradients [11, 12], and L-BFGS [13]. Though computationally tractable, the
quasi-Newtonian methods may suffer from convergence stagnation.

The Newton methods achieve faster and more reliable convergence, which catch
up with the nonlinearity through quadratic approximation of the target
function, see [14, 15, 16]. The often-used Gauss-Newton variant computes the
Hessian matrix approximately by linearizing the forward-problem operator. The
main disadvantage of the Newton methods applied to (1) is an enormous
computational burden stemming from the need to solve a large dense system of
linear equations at each Newton iteration. The system is solved iteratively,
usually with the conjugate gradients (CG). The matrix is ill-conditioned, and
therefore each Newton (external, nonlinear) iteration may require hundreds of
internal (linear) iterations. In addition, no efficient preconditioner is
known to improve this matrix’s conditioning. For these reasons, numerical
implementations of the Newton FWI are far less common than the quasi-Newton
ones.

It is natural to ask if we can obtain the fast Newton-like convergence with a
less resource-demanding procedure. The ultimate goal of this paper is to
demonstrate that the answer is yes.

Optimization (1) is known as the reduced-space approach in the optimization
community. Its alternative is the the full-space approach, also known as all-
at-once. In the reduced-space approach, only material parameters $\gamma$ are
optimized, whereas the connection between the wavefield in the medium and
$\gamma$ is prescribed by the forward-problem operator $\mathcal{F}$. In the
full-space approach, one simultaneously optimizes the control, state, and
adjoint-state variables. Upon formulating the necessary conditions and
discretizing the PDE, one obtains a nonlinear system of algebraic equations,
known as the Karush-Kuhn-Tucker (KKT) system, which is sparse. Its sparsity is
a clear advantage over formulation (1). The main challenge here is that the
problem becomes many times bigger than the reduced-space approach because the
state and adjoint-state variables for all sources are included in
optimization.

The full-space approach has been introduced to geophysics in the seminal paper
[17], followed by [18, 19]. Those authors cast the geophysical electromagnetic
inversion as a full-space optimization problem, derived the KKT system and
applied the Newton algorithm to it. Despite spectacular results, the full-
space approach had been considered impractical for years due to the enormous
size of the linear system, arising on each Newton step. Regarding FWI, the KKT
optimality conditions were derived in several papers [15, 16], but were viewed
as a theoretical tool rather than a basis for numerical implementation.

Significant progress has been made in this direction during the last decade.
Papers [20, 21, 22] report the application of the penalty method to the FWI,
marking an important step toward the full-space approach. The penalty method
has some theoretical drawbacks connected to the choice of the penalty
parameter [23] (but see recent developments [24]). The augmented Lagrangian
(AL) method is a viable alternative to the penalty method. A variant of it,
known as the alternating-direction method of multiplies (ADMM) [25], was
applied to the FWI with great success, see [26, 27]. In essence, the ADMM
tackles the full-space optimization by making alternating steps in the state,
adjoint-state, and control state spaces. Still, it is generally accepted that
the main weakness of the ADMM compared to the Newton iterations is its slow
convergence [25].

This paper studies the full-space optimization based on applying a Newton
solver to the original KKT matrix. The crucial component of our approach is a
suitable preconditioner based on a block-triangular factorization of the
permuted KKT matrix combined with approximations of individual blocks. Similar
ideas were exploited optimization with Maxwell’s equations [18] and the Stokes
equation [28]; thus, this paper extends this framework to Helmholtz’s
equation. We benchmarked our method against the standard reduced-space Gauss-
Newton FWI with the internal CG solver. In the numerical experiments conducted
the Marmousi model, we observed the 6$\times$ speedup. The method has the
potential for further acceleration. In particular, it will benefit from any
method designed to approximate the reduced-space Hessian, such as L-BFGS. The
central result of this paper is that it establishes the feasibility of Newton-
type optimization of the KKT system in application to seismic FWI.

The paper is organized as follows. In Section 2, we formulate the forward
problem. The reduced-space Gauss-Newton inversion is considered in Section 3.
Sections 4 and 5 are dedicated to the full-space approach and its iterative
solution. Numerical experiments are presented in Section 6. Concluding remarks
are given in Section 7.

##  2 The forward problem

In this section, we briefly review the forward problem. Let
$\Omega\subset\mathbb{R}^{2}$ be a bounded domain with boundary
$\partial\Omega$. Consider the 2D Helmholtz’s equation in the following form,

$\displaystyle-Z\frac{\partial}{\partial x}\left(\frac{1}{X}\frac{\partial
U}{\partial x}\right)-X\frac{\partial}{\partial
z}\left(\frac{1}{Z}\frac{\partial U}{\partial z}\right)-\omega^{2}\gamma
XZU=XZF,\text{ in }\Omega,$ (2) $\displaystyle U=0,\text{ on }\partial\Omega.$

Here $F(x,z)$ are the right-hand side, $U(x,z,\omega)$ are the acoustic field
(monochromatic pressure), $\omega$ is the circular frequency, $\gamma(x,z)$ is
the squared slowness, $\gamma=c^{-2}$, where $c(x,z)$ is the speed of sound.
Complex-valued functions of one variable $X(x)$ and $Z(z)$ implement the
perfectly-matched layers (PML). They deviate from 1.0 only in the damping
region on the sides of the computational domain to suppress the side
reflections, see Figure 1. On the PML implementation, see [29] and references
therein.

Figure 1: Schematic representation of the computational domain with a PML
region. In the PML region $X(x)\neq 1$, $Z(z)\neq 1$.

Let us rewrite (2) with a more compact notation as follows,

$\displaystyle\mathcal{D}U=F\text{ in }\Omega,$ (3) $\displaystyle U=0,\text{
on }\partial\Omega.$

where

$\mathcal{D}(\cdot):=-Z\frac{\partial}{\partial
x}\frac{1}{X}\frac{\partial}{\partial x}(\cdot)-X\frac{\partial}{\partial
z}\frac{1}{Z}\frac{\partial}{\partial z}(\cdot)-\omega^{2}\gamma XZ(\cdot)\,.$
(4)

Note that we replaced $XZF$ with $F$ in (3) because we assume that $F$ is
supported only in the core domain.

In a seismic survey, the acoustic field is excited by many monochromatic
sources independent of each other. Thus, we have a series of problems,

$\displaystyle\mathcal{D}U_{k}=F_{k},\text{ in }\Omega,$ (5) $\displaystyle
U_{k}=0,\text{ on }\partial\Omega,$

where index $k$ designates the source index. Let us combine solutions $U_{k}$
and right-hand sides $F_{k}$ into vector variables to make the formulas even
more compact. Starting from now on, $U$ is a direct product of all $U_{k}$,
$F$ is a direct product of $F_{k}$:

$U:=\prod_{k=0}^{K-1}U_{k},\quad F:=\prod_{k=0}^{K-1}F_{k},\quad k=0..K-1.$
(6)

With these notations, a set of $K$ problems (5) can be compactly written as

$\displaystyle\mathcal{D}U=F,\text{ in }\Omega,$ (7) $\displaystyle U=0,\text{
on }\partial\Omega.$

Now we introduce the observation operator $\mathcal{Q}$ which takes the
acoustic field and returns a vector of data $d$. We define $\mathcal{Q}$
through the convolution with the Dirac $\delta$-function (although any linear
functional can be used instead) as follows,

$d=\mathcal{Q}U:=\left\\{\int_{\Omega}\delta(x_{kj}-x)U_{k}\,dx\right\\},\quad
k=0,1,..K-1,\quad j=0,1,..N_{k}-1$ (8)

where $x_{kj}$ is the position of the $j$-th receiver of the $k$-th source,
$N_{k}$ is the number of receivers for the $k$-th source,
$d\in\mathbb{C}^{N}$, $N=\sum_{k=0}^{K-1}N_{k}$ is the total number of
receivers for all sources.

In what follows, we will use the discrete form of the forward problem. The
finite-difference (FD) discretization is assumed because it is used in our
numerical experiments. However, we might as well have performed discretization
by the finite-element of spectral-element methods.

Let $u_{k}$, $f_{k}$, and $s$ be finite-difference approximations of $U_{k}$,
$F_{k}$, and $\gamma$, respectively. We introduce the following notation,

$u=\prod_{k=0}^{K-1}u_{k},\quad f=\prod_{k=0}^{K-1}f_{k},\quad k=0..K-1.$ (9)

A set of discrete forward problems reads

$\begin{split}Au=f,\end{split}$ (10)

where $A$ is a block-diagonal matrix. Each diagonal block $A_{j}$ of matrix
$A$ is a sparse complex non-Hermitian matrix. If the PMLs are not used, then
$A_{j}$ turns into the real symmetric indefinite matrix,

$-\Delta-\omega^{2}M,$ (11)

where $\Delta$ is a discrete Laplacian, $M$ is a diagonal mass matrix,
containing values of $s$ on its diagonal $M=Is$ (it can be non-diagonal in
some FD schemes [30]).

##  3 The reduced-space Gauss-Newton inversion

In this section, we review the reduced-space Gauss-Newton inversion (RSGN).
Since we will eventually end up with the FD approximation, and thus discrete
variables will be grid functions, it is easier to formulate the inverse
problem in the discrete form (the discretize-then-optimize approach). In the
finite-element framework, a more convenient way would be to formulate the
optimality conditions in function spaces and then move on to discretization
(the optimize-then-discretize approach), see [31].

Let $S$ be a nonlinear operator that maps a given $s$ to $u$,

$u=S(s).$ (12)

Let $Q$ be a finite-dimensional operator that maps a given discrete acoustic
field $u$ to data $d$,

$d=Qu.$ (13)

Discrete coefficient $s$ is iterated as $s_{n+1}=s_{n}+\,\delta s\,$, where
$\,\delta s\,$, known as the model update. On each iteration, $s$ is a
solution to the following minimization problem,

$\underset{\delta s}{\text{minimize}}\quad\frac{1}{2}\left\lVert
W(r-QG\,\delta s\,)\right\rVert^{2}_{2}+\frac{\varepsilon}{2}\left\lVert
s_{n}+\,\delta s\,\right\rVert_{2}^{2}.$ (14)

Here $W$ is a real-valued diagonal matrix of weights designed primarily to
balance signal attenuation with offset from the source. Vector $r$ is the data
residual,

$r=d-QS(s_{n}),$ (15)

$G$ is a matrix of the Frechét operator of $S$ at point $s=s_{n}$, e.i.

$G=\left.\frac{\partial S}{\partial s}\right|_{s=s_{n}}.$ (16)

Solution to (14) satisfies

$\left(G^{*}Q^{*}W^{T}WQG+L\right)\,\delta s\,=W^{T}G^{*}Q^{*}r-Ls_{n},$ (17)

where $L=\varepsilon I$. The following formula defines matrix $G$,

$G=A^{-1}P,$ (18)

where $P$ is a diagonal matrix containing forward solution $u_{n}$ at
$s=s_{n}$,

$P:=\omega^{2}Iu_{n}.$ (19)

To verify (18),(19), one can subtract the two discrete Helmholtz’s equations

$\displaystyle\left(-\Delta-\omega^{2}I\,(s_{n}+\,\delta
s\,)\right)(u_{n}+\delta u)=f,$ (20)
$\displaystyle\left(-\Delta-\omega^{2}I\,s_{n}\right)u_{n}=f,$

although the derivation can be conducted directly in the Hilbert spaces. In
geophysics, the matrix

$J:=QA^{-1}P,$ (21)

is known as the Jacobian, the matrix

$H:=J^{*}W^{T}WJ+L$ (22)

is the (regularized) Hessian, and vector

$g:=W^{T}J^{*}r-Ls_{n}$ (23)

is called the gradient. Using this definitions, we rewrite (17) as

$H\,\delta s\,=g.$ (24)

For any $\varepsilon>0$ matrix $H$ is positive-definite, and thus (24) has a
unique solution. It should be noted that, even if $\varepsilon=0$ (singular
$H$), then a Krylov solver will deliver a meaningful estimate of $\,\delta
s\,$, that is, a projection to a corresponding Krylov subspace.

The system matrix in (24) is dense and large. In real-life 3D problems, the
system (24) can be solved only iteratively, with the CG being a natural
choice. We will abbreviate this algorithm as RSGN-CG. The CG multiplies the
system matrix $J^{*}J$ by a vector once per iteration. Because of (21) it
implies two linear solves per CG iteration - the forward simulation with $A$
and the adjoint simulation with $A^{*}$. We can rewrite (24) as

$\left(\sum_{j=0}^{K-1}J^{*}_{j}W_{j}^{T}W_{j}J_{j}+L\right)\,\delta
s\,=\sum_{j=0}^{K-1}W_{j}^{T}J^{*}_{j}r_{j}-Ls_{n}.$ (25)

The forward and adjoint simulations for different sources can be performed in
parallel but these results must be combined at each CG iteration.

However, the most critical problem is a very slow convergence of the CG. The
system matrix (24) has a very high condition number for all meaningful values
of $\varepsilon$. Thus, tens or hundreds of forward simulations (for all
sources) may be needed to generate a good update $\,\delta s\,$. In addition,
preconditioning of a regularized normal system of linear equitations is
notoriously difficult.

##  4 The full-space Gauss-Newton inversion

In this section, we describe the full-space Gauss-Newton (FSGN) inversion. We
start by noting that the discrete optimization problem (14) is equivalent to
the following constrained optimization problem,

$\begin{split}\underset{\,\delta
s\,}{\text{minimize}}\quad\Psi:=\frac{1}{2}\left\lVert W\left(r-Q\,\delta
u\,\right)\right\rVert^{2}_{2}+\frac{\varepsilon}{2}\left\lVert s_{n}+\,\delta
s\,\right\rVert_{2}^{2},\\\ \text{ s.t. }\quad A\,\delta u\,=P\,\,\delta
s\,.\end{split}$ (26)

Here $n$ is the Newton iteration number, $\,\delta u\,$ is the update to
$u_{n}$, e.i. $u_{n+1}\approx u_{n}+\,\delta u\,$. The target functional
$\Psi$ is quadratic and bounded away from zero for any $\varepsilon>0$. The
constraints are linear. For such problems, known as linear-quadratic, the
existence of unique state and control variables are readily established, for
example, [32][Thm 1.43]. Moreover, since we start from the discrete
formulation, the solvability is proved simply by establishing the equivalence
between the matrix form and the normal system (24), see below.

Let us form the Lagrangian,

$\mathcal{L}:=\Psi+\lambda^{*}(A\,\delta u\,-P\,\,\delta s\,),$ (27)

where $\lambda$ is the adjoint-state variable (a complex-valued grid
function). By applying the optimality conditions [32], we obtain the KKT
system,

$\begin{split}-Q^{*}W^{T}W(r-Q\,\delta u\,)+A^{*}\lambda=0,\\\
L(s_{n}+\,\delta s\,)-P^{*}\lambda=0,\\\ A\,\delta u\,-P\,\delta
s\,=0.\end{split}$ (28)

This system can be rewritten in the matrix form,

$\underbrace{\left[\begin{array}[]{ccc}F&O&A^{*}\\\ O&L&-P^{*}\\\ A&-P&O\\\
\end{array}\right]}_{\mathcal{M}}\underbrace{\left[\begin{array}[]{c}\,\delta
u\\\ \,\delta s\\\
\lambda\end{array}\right]}_{\xi}=\underbrace{\left[\begin{array}[]{c}Q^{*}W^{T}Wr\\\
-Ls_{n}\\\ 0\end{array}\right],}_{b}$ (29)

where $F=Q^{*}W^{T}WQ$.

Lets us prove that (29) is equivalent to (24) for it establishes solvability
of the KKT system. Using the third equation of (29), we compute $\,\delta
u\,$,

$\,\delta u\,=A^{-1}P\,\delta s\,.$ (30)

Using this result, we compute $\lambda$ from the first equation,

$\lambda=A^{-*}\left(Q^{*}r-Q^{*}QA^{-1}P\,\delta s\,\right),$ (31)

where notation $A^{-*}$ means $(A^{*})^{-1}$. Now we substitute (30) and (31)
into the second equation of (29) and get

$\left(P^{*}A^{-*}Q^{*}QA^{-1}P+L\right)\,\delta
s\,=P^{*}A^{-*}Q^{*}r-Ls_{n}.$ (32)

From (21) we see that (32) is equivalent to (24).

From the preceding discussion, it is clear that the model update $\,\delta
s\,$ computed from either the KKT system (29) or the normal system (24) is
identical, provided that the linear systems are solved exactly. The difference
between the two approaches lies in the structure of the system matrix. The
system matrix in (24) is complex-valued, dense, Hermitian positive-definite,
and relatively small. For small- to medium-size 3D problems with rigorously
selected FD grids and frequencies, satisfactory results have been obtained
with the use of the sparse LU decomposition of the system $A$ [8, 33]. In
general, however, a parallel on-the-fly application of $A^{-1}$ and $A^{-*}$
is the only option. The system matrix in (29) is sparse complex-valued
Hermitian indefinite. It is many times larger than the normal matrix, and it
precludes from applying direct solvers in any form. Its iterative solution is
considered in the next section.

##  5 Iterative solution of the KKT system

In this section, we design an iterative solver for the KKT system (29). The
system matrix $\mathcal{M}$ is can be explicitly partitioned across sources as
follows,

$\left[\begin{array}[]{@{}c|c|c@{}}\begin{matrix}F_{0}&&\\\ &\ddots&\\\
&&F_{K-1}\\\ \end{matrix}&O&\begin{matrix}A^{*}_{0}&&\\\ &\ddots&\\\
&&A^{*}_{K-1}\\\ \end{matrix}\\\ \hline\cr
O&L&\begin{matrix}-P^{*}_{0}\ldots-P^{*}_{K-1}\end{matrix}\\\
\hline\cr\begin{matrix}A_{0}&&\\\ &\ddots&\\\ &&A_{K-1}\\\
\end{matrix}&\begin{matrix}-P_{0}\\\ \vdots\\\
-P_{K-1}\end{matrix}&O\end{array}\right]\left[\begin{array}[]{c}\,\delta
u\,_{0}\\\ \vdots\\\ \,\delta u\,_{K-1}\\\ \hline\cr\,\delta s\\\
\hline\cr\lambda_{0}\\\ \vdots\\\
\lambda_{K-1}\end{array}\right]=\left[\begin{array}[]{c}Q_{0}^{*}W_{0}^{T}W_{0}r_{0}\\\
\vdots\\\ Q_{K-1}^{*}W_{K-1}^{T}W_{K-1}r_{K-1}\\\ \hline\cr-Ls_{n}\\\
\hline\cr 0\\\ \vdots\\\ 0\end{array}\right].$ (33)

Block $F$ is a positive semidefinite matrix, which is very sparse. If
receivers are located at the grid nodes, then $F$ is a diagonal matrix where
the number of non-zero entries in each $F_{j}$ block equals the number of data
points for the $j$-th source. Block $A$ is a block-diagonal matrix of full
rank whose spectrum lies in the left and right half-planes. Block $L$ is
positive-definite here but may become positive semidefinite if $L$ represents
the Laplacian. In our case, $L=\varepsilon I$ and so its diagonal. Block $P$
is a rectangular matrix consisting of sparse blocks. All submatrices of $P$
are diagonal; see definition (19). The pattern of the matrix $\mathcal{M}$ for
a test problem consisting of 25 sources is presented in Figure 2.

Figure 2: Partitioning (red lines) and non-zero elements (black dots) of the
KKT matrix for a test problem. The problem contains 25 sources, each system
matrix is 120,000 $\times$ 120,000. The size of the $L$ block is 72,000
$\times$ 72,000.

Our preconditioner is based on a block-triangular decomposition, which is one
of the most efficient strategies for preconditioning saddle-point problems. In
the context of the KKT matrix, this approach was considered in [18, 28].

The matrix $\mathcal{M}$ has singular (1,1) and (3,3) blocks, and so the
direct block-LU decomposition is not possible. Thus, we perform block-LU
factorization of the permuted matrix as follows,

$\left[\begin{array}[]{ccc}A&O&-P\\\ F&A^{*}&O\\\ O&-P^{*}&L\\\
\end{array}\right]=\left[\begin{array}[]{ccc}I&O&O\\\ FA^{-1}&I&O\\\
O&-P^{*}A^{-*}&I\\\ \end{array}\right]\left[\begin{array}[]{ccc}A&O&-P\\\
O&A^{*}&FA^{-1}P\\\ O&O&H\\\ \end{array}\right],$ (34)

where $H=L+P^{*}A^{-*}FA^{-1}P$. The $U$-factor is an excellent preconditioner
since applying its inverse to the initial (permuted) matrix results in a
matrix whose spectrum consists of a single eigenvalue equal to 1. This
approach would be equivalent to the reduced-space method and has no advantages
because it requires a linear solve with the reduced Hessian $H$.

Let us approximate (34) by assuming $F\approx 0$, which implies $H\approx L$.
Plugging $F=0,H=L$ into (34), we get

$\left[\begin{array}[]{ccc}I&O&O\\\ O&I&O\\\ O&-P^{*}A^{-*}&I\\\
\end{array}\right]\left[\begin{array}[]{ccc}A&O&-P\\\ O&A^{*}&O\\\ O&O&L\\\
\end{array}\right]=\left[\begin{array}[]{ccc}A&O&-P\\\ O&A^{*}&O\\\
O&-P^{*}&L\\\ \end{array}\right].$ (35)

By reversing the permutation, we define the preconditioner to the original
system matrix $\mathcal{M}$ as follows,

$\mathcal{P}:=\left[\begin{array}[]{ccc}O&O&A^{*}\\\ O&L&-P^{*}\\\ A&-P&O\\\
\end{array}\right].$ (36)

The preconditioned matrix , $\mathcal{P}^{-1}\mathcal{M}$ is non-symmetric;
therefore, a non-symmetric Krylov solver is required to solve it. We selected
the GMRes and abbreviated the resulting method as FSGN-GMRes.

The preconditioner (36) has several appealing properties. First, it requires
three forward solves with matrices $A$, $A^{*}$, and $L$, e.i. almost the same
work as needed for a single iteration of the CG with a reduced Hessian.
Second, the $A$ and $A^{*}$ both split into $K$ independent linear systems.
Thus, the solutions with blocks (1,3) and (3,1) can be performed in parallel.
Third, the preconditioner performs approximation $H\approx L$, but other
approximations to $H$ are possible. In particular, one can use the L-BFGS
approximation in place of block $L$, see [28].

Now, we define another preconditioner, which is obtained from $\mathcal{P}$ by
approximating $A$ by another matrix $\tilde{A}$ as follows,

$\tilde{\mathcal{P}}:=\left[\begin{array}[]{ccc}O&O&\tilde{A}^{*}\\\
O&L&-P^{*}\\\ \tilde{A}&-P&O\\\ \end{array}\right].$ (37)

In this paper, we utilized the incomplete LU (ILU) factorization to
approximate $A$. This choice is motivated by the fact that the ILU is widely
available and requires no tuning apart from selecting the fill-in level.

We will demonstrate below that the FSGN-GMRes with the $\tilde{\mathcal{P}}$
preconditioner is superior to the RSGN-CG. Still, it is important to recognize
that more efficient strategies to approximate $A$ are likely to exist. In
particular, the shifted-Laplacian preconditioner leveraged with the multigrid
[34] deserves a dedicated study.

##  6 Numerical experiments

This section is dedicated to numerical experiments. The optimization part of
FWI was programmed in Python, whereas the forward and adjoint simulation
problem is coded in C++. The forward problem was solved directly with a sparse
LU factorization. All calculations were serial.

We conducted numerical experiments using a part of the Marmousi velocity model
(Figure 3). This P-wave velocity distribution was used in the data simulation,
representing the geologic medium we want to reconstruct.

Figure 3: The distribution of the P-wave velocity, $c$, used to simulate
synthetic data. The numbers on the color scale are the velocity in km/s. The
maroon lines on the sides indicate the boundaries of the PML region.

The purpose of the first experiment was to compare, inside a single GN step,
the performance of the CG applied to the normal system matrix with the
preconditioned GMRes applied to the KKT matrix. With the GMRes, we checked two
preconditioners, $\mathcal{P}$ and $\tilde{\mathcal{P}}$. The velocity
distribution shown in Figure 4(a) was the initial model for the inversion. The
computed slowness updates are compared in Figure 4(b)-(d).

---
(a)
(b)
(c)
(d)
Figure 4:  Comparison of the the two slowness updates computed at the end of a
single GN iteration. Panel (a) shows the initial velocity distribution. Panel
(b) depicts the slowness update obtained by applying the CG solver to the
normal system of linear equations. Panel (c) shows the slowness update by the
GMRes solver with $\mathcal{P}$ applied to the KKT system. Panel (d) same as
(c) but with preconditioner $\tilde{\mathcal{P}}$ Note that in (a) color
represents the P-wave velocity in km/s, whereas in (b)-(d) color represents
update to the squared slowness in $s^{2}/m^{2}$.

The preconditioner $\tilde{\mathcal{P}}$ was created by making use of ILU with
the fill-in parameter 21. It is not not possible to compare the CG and GMRes
directly, because the CG computes $E_{cg}:=\left\lVert H\,\delta
s\,_{n}-g\right\rVert_{2}$ at each iteration, whereas GMRes computes
$E_{gmres}:=\left\lVert\mathcal{M}\xi_{n}-b\right\rVert_{2}$. To make the
comparison, we separately computed values of $E_{cg}$ from the GMres iterates
$\xi_{n}$. The convergence of the iterative methods are presented in Figure 5.

Figure 5: Evolution of data misfit for the CG applied to the normal system and
the preconditioned GMRes applied to the KKT matrix. The $Y$-axis represents
the quantity $\left\lVert H\,\delta s\,_{n}-g\right\rVert_{2}/\left\lVert
g\right\rVert_{2}$.

The CG and GMRes achieved approximately the same tolerance on iterations 19
and 30, respectively (Figure 5). As expected, all methods generated almost
identical updates. The convergence rate of the CG is the same as that of the
GMRes with the $\mathcal{P}$, although the GMRes happened to reach the desired
tolerance by 3 iterations faster. The convergence of the GMRes with
$\tilde{\mathcal{P}}$ is slower than that with $\mathcal{P}$, but the cost is
much lesser. The performance of the algorithms is compared in Table 1. We see
the dramatic acceleration of matrix-vector multiplication in the GMRes
compared to the CG. This is because the linear solves with factorized matrices
$\tilde{A}$ and $\tilde{A}^{*}$ inside $\tilde{\mathcal{P}}$ is much cheaper
than the linear solves with $A$ and $A^{*}$ during the CG matrix-vector
multiplication. Another observation is that we needed a high value of fill-in
(we used 21) to have an acceptable approximation $A\approx\tilde{A}$.

Table 1: Linear solvers comparison. The column ’CG’ refers to the the CG applied to the normal system of linear equations, ’GMRes-full’ is the GMRes applied to the KKT matrix preconditioned with $\mathcal{P}$, ’GMRes-approx’ is the GMRes preconditioned with $\tilde{\mathcal{P}}$. | CG | GMRes-full | GMRes-approx
---|---|---|---
Iterations | 19 | 16 | 30
Time per iteration, s | 264 | 267 | 18
ILU initialization, s | 0 | 0 | 121
Total time, s | 4800 | 5391 | 834

In the second numerical experiment, we demonstrate that our approach serves as
an efficient engine for the multi-frequency FWI. The data were simulated at
frequencies from 5 to 40 Hz with a step of 2.5 Hz. At each frequency in the
sequence, from lower to higher, a single GN step was performed, and the
updated model was passed to the next frequency. The regularization parameter
$\varepsilon$ was changing linearly from 10 to 1E5. The number of internal
GMRes iterations was set to 30. As indicated in Figure 6, the data misfit was
improved after the inversion. The improvement becomes more pronounced at
higher frequencies because the data magnitudes increase with the frequency and
the initial model describes high-frequency data very poorly.

Figure 6: The relative residual misfit on Newton iterations: ’resid-norm-ini’
is $\left\lVert d_{initial}-d\right\rVert_{2}/\left\lVert d\right\rVert_{2}$,
’resid-norm-fin’ is $\left\lVert d_{final}-d\right\rVert_{2}/\left\lVert
d\right\rVert_{2}$.

The initial velocity model, several intermediate models, and the final model
are presented in Figure 7. In this numerical experiment, we set the term
$-Ls_{n}$ in (29) to zero and, thus, the algorithm was penalizing the norm of
the update $\delta s$, not the updated model $s_{n}+\delta s$. This was the
reason for the high-velocity build-up at $X=2500$. We stress again that
identical FWI results would have been obtained if we had applied the standard
RSGN-CG optimization (provided, of course, that the values of error $E_{cg}$
match in both approaches). The difference between the two approaches is in
their performance.

---
(a) Initial
(b) 5 Hz
(c) 10 Hz
(d) 17.5 Hz
(e) 25 Hz
(f) 40 Hz
Figure 7: The initial velocity distribution (a) and updated distributions
(b)-(f) at different frequencies. The color scales here and in Figure 3 are
identical. The X- and Z-axis are marked in meters. Similar, if not identical,
updates would have been obtained if we had performed the standard RSGN-CG
optimization with the same accuracy of data fit.

##  7 Conclusions

This paper proposes a new numerical method for the seismic FWI, based on the
full-space approach. Prior work has demonstrated that the full-space approach
can be used as a framework for efficient waveform inversion algorithms. Those
works have focused either on the penalty method with the primal-dual descend
optimization or the AL approach with the ADMM optimization procedure. These
approaches have demonstrated impressive results, but may suffer from the slow
convergence compared to the Newton iterations [25].

This paper studies the iterative full-space FWI by the Newton method. A
preconditioned Krylov solver is applied to the KKT system on each Newton
iteration. Our method computes the material parameter updates identical to
those calculated by inverting the reduced Hessian, provided that the linear
systems are solved exactly. The critical element of our approach is a special
preconditioner. The preconditioned iterative algorithm requires two forward
simulations (exact or approximate) at each iteration plus an inexpensive solve
with the constraint matrix $L$. We considered two variants of the
preconditioner: preconditioner $\mathcal{P}$ computes linear solves with $A$
and $A^{*}$ exactly, whereas preconditioner $\tilde{\mathcal{P}}$ computes the
inexact solves. A single application of $\mathcal{P}$ should therefore have
almost the same cost compared to a single iteration of the CG applied to the
normal system. Our numerical experiments indicated that these two methods
indeed have very close performance. On the other hand, a single application of
$\tilde{\mathcal{P}}$ is much cheaper. In our numerical experiments, this
approach reduced the execution time considerably.

For a real-scale 3D FWI, several problems need to be addressed. First, the
forward and adjoint simulations need to be based on the theory of linear
elasticity, so it remains to be shown that our approach works for elastic FWI.
Second, direct linear solvers become unfeasible in the 3D applications, even
for the acoustic equations, not to mention the linear elasticity. Thus, a
reliable preconditioned iterative solver should be applied for the forward and
adjoint simulations. Substantial progress in this direction has been made; see
[29] and references therein. In principle, any preconditioner designed for the
forward simulation can be used inside $\tilde{\mathcal{P}}$. Third, the
convergence can be improved by a more sophisticated approximation of the
reduced Hessian, such as the L-BFGS or the background Hessian. Virtually any
method designed to accelerate the reduced-space Newton method can be used to
compute the $H$-block solve within our approach. Finally, using the BiCGStab
instead of the GMRes will eliminate the memory needed to store GMRes iterates
between restarts. Future work, therefore, should include a study of these
options.

##  8 Acknowledgments

This research was partially supported by the Russian Science Foundation,
project no. 21-11-00139. The authors acknowledge computational resources
granted by Complex for Simulation and Data Processing for Mega-Science
Facilities at NRC ”Kurchatov Institute”, http://ckp.nrcki.ru/.

## References

  * [1] P. Lailly.  The seismic inverse problem as a sequence of before stack migrations.  In Conference on Inverse Scattering, Theory and Application, pages 206–220. SIAM, 1983.
  * [2] A. Tarantola.  Inversion of seismic reflection data in the acoustic approximation.  Geophysics, 49(8):1259–1266, 1983.
  * [3] P. R. Mora.  Nonlinear two-dimensional elastic inversion of multi-offset seismic data.  Geophysics, 52(9):1211–1228, 1987.
  * [4] J. Virieux and S. Operto.  An overview of full-waveform inversion in exploration geophysics.  Geophysics, 74:WCC1–WCC26, 2009.
  * [5] J. Tromp.  Seismic wavefield imaging of Earth’s interior across scales.  Nature Reviews Earth & Environment, 1:40–53, 2020.
  * [6] R. Plessix.  Three-dimensional frequency-domain full-waveform inversion with an iterative solver.  Geophysics, 74(6):WCC149–WCC157, 2009.
  * [7] M. Warner, A. Ratcliffe, T. Nangoo, J. Morgan, A. Umpleby, N. Shah, V. Vinje, I. Štekl, L. Guasch, C. Win, G. Conroy, and A. Bertrand.  Anisotropic 3d full-waveform inversion.  Geophysics, 78(2):R59–R80, 2013.
  * [8] S. Operto, A. Miniussi, R. Brossier, L. Combe, L. Métivier, V. Monteiller, A. Ribodetti, and J. Virieux.  Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation.  Geophysical Journal International, 202(2):1362–1391, 07 2015.
  * [9] S Operto and A Miniussi.  On the role of density and attenuation in three-dimensional multiparameter viscoacoustic VTI frequency-domain FWI: an OBC case study from the North Sea.  Geophysical Journal International, 213(3):2037–2059, 2018.
  * [10] A. Fichtner, B. L. N. Kennett, H. Igel, and H.-P. Bunge.  Full seismic waveform tomography for upper-mantle structure in the australasian region using adjoint methods.  Geophysical Journal International, 179(3):1703–1725, 2009.
  * [11] C. Shin, S. Jang, and D.-J. Min.  Improved amplitude preservation for prestack depth migration by inverse scattering theory.  Geophysical Prospecting, 49(5):592–606, 2001.
  * [12] W.A. Mulder and R.-E. Plessix.  Exploring some issues in acoustic full waveform inversion.  Geophysical Prospecting, 56(6):827–841, 2008.
  * [13] R. Brossier, S. Operto, and J. Virieux.  Seismic imaging of complex onshore structures by 2d elastic frequency-domain full-waveform inversion.  Geophysics, 74(6):WCC105–WCC118, 2009.
  * [14] R. G. Pratt, C. Shin, and G. J. Hick.  Gauss–Newton and full Newton methods in frequency–space seismic waveform inversion.  Geophysical Journal International, 133(2):341–362, 1998.
  * [15] I. Epanomeritakis, V. Akçelik, O Ghattas, and J Bielak.  A Newton-CG method for large-scale three-dimensional elastic full-waveform seismic inversion.  Inverse Problems, 24(3):034015, 2008.
  * [16] L. Métivier, R. Brossier, S. Operto, , and J. Virieux.  Full waveform inversion and the truncated newton method.  SIAM Review, 59(1):153–195, 2017.
  * [17] E.Haber, U. M. Ascher, and D. Oldenburg.  On optimization techniques for solving nonlinear inverse problems.  Inverse Problems, 16(4):1263, 2000.
  * [18] E. Haber and U. M. Ascher.  Preconditioned all-at-once methods for large, sparse parameter estimation problems.  Inverse Problems, 17:1847–1864, 2001.
  * [19] E. Haber, U. M. Ascher, and D. W. Oldenburg.  Inversion of 3D electromagnetic data in frequency and time domain using an inexact all-at-once approach.  Geophysics, 69(5):1216–1228, 2004.
  * [20] A. Abubakar, W. Hu, T. M. Habashy, and P. M. van den Berg.  Application of the finite-difference contrast-source inversion algorithm to seismic full-waveform data.  Geophysics, 74(6):WCC47–WCC58, 2009.
  * [21] T. van Leeuwen and F. J. Herrmann.  Mitigating local minima in full-waveform inversion by expanding the search space.  Geophysical Journal International, 195(1):661–667, 2013.
  * [22] T. van Leeuwen and F. J. Herrmann.  A penalty method for PDE-constrained optimization in inverse problems.  Inverse Problems, 32:015007, 2016.
  * [23] J. Nocedal and S. J. Wright.  Numerical Optimization.  Springer, New York, NY, USA, 2e edition, 2006.
  * [24] G. Rizzuti, M. Louboutin, R. Wang, and F. J. Herrmann.  A dual formulation of wavefield reconstruction inversion for large-scale seismic inversion.  Geophysics, 86(6):R879–R893, 2021.
  * [25] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Found. Trends Mach. Learn., 3(1):1–122, jan 2011.
  * [26] H. S. Aghamiry, A. Gholami, and S. Operto.  Improving full-waveform inversion by wavefield reconstruction with the alternating direction method of multipliers.  Geophysics, 84(1):R125–R148, 2019.
  * [27] K. Aghazade, A. Gholami, H. S. Aghamiry, and S. Operto.  Anderson-accelerated augmented lagrangian for extended waveform inversion.  Geophysics, 87(1):R79–R91, 2022.
  * [28] G. Biros and O. Ghattas.  Parallel Lagrange–Newton–Krylov–Schur methods for pde-constrained optimization. part i: The krylov–schur solver.  SIAM Journal on Scientific Computing, 27(2):687–713, 2005.
  * [29] N. Yavich, N. Khokhlov, M. Malovichko, and M.S Zhdanov.  Contraction operator transformation for the complex heterogeneous helmholtz equation.  Computers and Mathematics with Applications, 86:63–72, 2021.
  * [30] B. Hustedt, S. Operto, and J. Virieux.  Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modelling.  Geophysical Journal International, 157(3):1269–1296, 2004.
  * [31] M.S. Malovichko, A.V. Tarasov, N.B. Yavich, and K.V. Titov.  Application of optimal control to inversion of self-potential data: theory and synthetic examples.  IEEE Transactions on Geoscience and Remote Sensing, in press:1–1, 2021.
  * [32] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich.  Optimization with PDE Constraints.  Springer, 2009.
  * [33] V. Kostin, S. Solovyev, A. Bakulin, and M. Dmitriev.  Direct frequency-domain 3D acoustic solver with intermediate data compression benchmarked against time-domain modeling for full-waveform inversion applications.  Geophysics, 84(4):T207–T219, 2019.
  * [34] Y.A. Erlangga, C. Vuik, and C.W. Oosterlee.  On a class of preconditioners for the helmholtz equation.  Applied numerical mathematics, (50):409–425, 2004.