#  Tele-Correlation: Calibrating Shear-Shear Correlation with Real Data

Zhi Shen  Jun Zhang  Cong Liu  Hekun Li  Haoran Wang  Zhenjie Liu  Jiarui Sun

###### Abstract

Tele-correlation refers to the correlation of galaxy shapes with large angular
separations (e.g., $>100$ degrees). Since there are no astrophysical reasons
causing such a correlation on cosmological scales, any detected tele-
correlation could disclose systematic effects in shear-shear correlation
measurement. If the shear estimators are measured on single exposures, we show
that the field distortion (FD) signal associated with the galaxy position on
the CCD can be retained and used in tele-correlation to help us directly
calibrate the multiplicative and additive biases in shear-shear correlations.
We use the DECaLS shear catalog produced by the Fourier_Quad pipeline to
demonstrate this idea. To our surprise, we find that significant
multiplicative biases can arise (up to more than 10%) due to redshift binning
of the galaxies. Correction for this bias leads to about 1$\sigma$ increase of
the best-fit value of $S_{8}$ from $0.760^{+0.015}_{-0.017}$ to
$0.777^{+0.016}_{-0.019}$ in our tomography study.

##  1 Introduction

Weak lensing has been established as a powerful probe of the cosmic structure.
Many large scale galaxy surveys set weak lensing as their primary scientific
goals, including Canada-France-Hawaii Telescope Lensing Survey
(CFHTLenS111www.cfhtlens.org/)[1], Dark Energy Survey
(DES222www.darkenergysurvey.org) [2], Hyper Suprime-Cam Subaru Strategic
Program (HSC333hsc.mtk.nao.ac.jp/ssp/)[3], Kilo-Degree Survey
(KiDS444kids.strw.leidenuniv.nl)[4]. In these studies, it has been found that
the parameter $S_{8}$ from weak lensing is smaller than that from the cosmic
microwave background measurement [5] at about $2-3\sigma$ level. To firmly
establish this result, however, we need to carefully examine all possible
sources of systematic errors, which is well known to be difficult. We expect
these issues to be better addressed in Stage IV surveys, including Euclid [6],
the Large Synoptic Survey Telescope555www.lsst.org/(LSST), the China Space
Station Telescope (CSST) [7], and Roman [8, 9], all of which are going to
observe billions of galaxy images for accurate weak lensing measurement.

So far, major weak lensing surveys calibrate the systematics using image
simulations. For example, CFHTLenS takes image simulations to quantify the
dependence of shear bias on the signal-to-noise-ratio (SNR) and galaxy size
[10, 11]; HSC calibrates the multiplicative and additive bias as a function of
SNR and image resolution in simulations [12]; KiDS team calibrates the shear
catalog with a comprehensive model considering the SNR, the galaxy size and
ellipticity, as well as the size and ellipticity of the point spread function
(PSF) [13]; The DES team corrects the multiplicative bias with a multiple
effective redshift distribution, which is achieved using simulations [14, 15].

Ideally, one would hope to test the shear recovery accuracy directly on the
real data. One way to do so is to make use of the field distortion (FD) signal
that naturally exists in all optical images, as proposed in [16]. By grouping
the galaxies according to their underlying FD signals, and observing how well
the galaxy shear estimators can recover the FD (measured in astrometry), one
can get a direct estimate of the multiplicative and additive shear biases.
This idea has been successfully applied to the processing of the CFHTLenS data
[16] and the DECaLS data [17] with the Fourier_Quad shear measurement method,
and has helped discover a selection effect due to the presence of image
boundaries [18]. This kind of test automatically involves all observational
and instrumental effects, and has been proved as a robust way for calibrating
shape measurement.

In this work, we show that the FD test can be extended to calibrate potential
biases for the two-point statistics, i.e., shear-shear correlations, which is
commonly believed to be more challenging to measure due to its small amplitude
($\lesssim 10^{-4}$). We note that this is not completely equivalent to the
one-point test, because even if the underlying shear error is close to zero on
average, it may still have non-zero spatial correlations due to residual
systematic effects on the focal plane.

The idea of the two-point (2pt hereafter) FD test is to cross-correlate the
shear estimators of two remotely separated galaxies, which we call Tele-
Correlation (TC hereafter). Such a correlation should not contain
contributions from the astrophysical signals. If the galaxy shear estimator is
measured on each exposure individually, so that the FD shear information is
kept in the catalog for each galaxy image, the galaxy pairs involved in the TC
measurement can be grouped according to the products of their underlying FD
shears, forming a straightforward way of calibrating for the potential biases
in the shear-shear correlation.

In this paper, we apply the Tele-Correlation method on the DECaLS shear
catalog produced by the Fourier_Quad (FQ hereafter) shear measurement method
[17]. In §2 we introduce the DECaLS shear catalog, and the methods we use for
the shear estimation and shear-shear correlation. In §3 we introduce the
concept of TC, and show the results from TC, mainly including the shear biases
from each redshift bin for our tomographic study. We show how the biases
change our cosmology constraints in §4. Finally, we give our main conclusion
in §5.

##  2 Data and Shear Catalog

Our shear catalog is based on the imaging data of the Dark Energy Camera
Legacy Survey (DECaLS)666www.legacysurvey.org. The total sky coverage is about
10,000 $\mathrm{deg}^{2}$, taken by the Dark Energy Camera (DECam) on the
Blanco 4 m telescope of the Cerro Tololo Inter-American Observatory. The image
files are pre-processed through the “Community Pipeline” to remove the
instrumental effects, and the sky backgrounds are kept. The shear catalog is
obtained using the Fourier_Quad pipeline, which evolves from the FQ shear
measurement method [19, 20, 21, 22, 23, 16, 17]. It involves all the necessary
steps for achieving accurate galaxy shape measurement, including background
removal, astrometric calibration, PSF reconstruction, etc.. The pipeline is
applied to the data of all three bands: g, r, z. With the field distortion
test, it is found that the quality of the z-band shear catalog is much better
than those of the other two bands. We therefore only use the z-band data in
this work for the shear-shear correlation. The typical galaxy number density
is about 3–5 per square arcmin [17]. The photo-$z$ we use in this paper is
measured by [24], and for cross-checking the redshift, we also use the
photo-$z$ catalog from [25]. In the following tests, we set a sample cut with
those criteria: we choose those galaxies from the z-band with signal-to-noise
ratio $\nu_{F}$ (defined in Fourier space, see [26]) is larger than 4.
Additionally, we exclude galaxies with large field distortion signals:
$|g^{f}_{1,2}|>0.0015$, which are located near the edge of the exposure. We
also exclude galaxies obtained from the problematic chips and those with the
z-band magnitude larger than 21. After these cuts, there are about one galaxy
per square arcmin. Given that the FQ shear estimators are measured on
individual exposures, we only use the shear estimators from two different
exposures to measure their correlations. This is for the purpose of avoiding
possible systematic errors from correlated PSF residuals [27].

The FQ shear measurement method uses the multipole moments of the galaxy power
spectrum to form the shear estimators. They are defined as:

$\displaystyle G_{1}$ $\displaystyle=-\frac{1}{2}\int
d^{2}\vec{k}(k_{x}^{2}-k_{y}^{2})T(\vec{k})M(\vec{k})$ (2.1) $\displaystyle
G_{2}$ $\displaystyle=-\int d^{2}\vec{k}k_{x}k_{y}T(\vec{k})M(\vec{k})$
$\displaystyle N$ $\displaystyle=\int
d^{2}\vec{k}[k^{2}-\frac{\beta^{2}}{2}k^{4}]T(\vec{k})M(\vec{k})$
$\displaystyle U$ $\displaystyle=-\frac{\beta^{2}}{2}\int
d^{2}k(k_{x}^{4}-6k_{x}^{2}k_{y}^{2}+k_{y}^{4})T(\vec{k})M(\vec{k})$
$\displaystyle V$ $\displaystyle=-2\beta^{2}\int
d^{2}k(k_{x}^{3}k_{y}-k_{x}k_{y}^{3})T(\vec{k})M(\vec{k}),$

in which $M(\vec{k})$ is the 2D galaxy power spectrum corrected by terms
related to the background noise and the Poisson noise [22]. $T(\vec{k})$ is
the factor for converting the PSF to a Gaussian form, i.e.:

$T(\vec{k})={\big{|}\widetilde{W}_{\beta}(\vec{k})\big{|}}^{2}/{\big{|}\widetilde{W}_{PSF}(\vec{k})\big{|}}^{2}$
(2.2)

in which $\widetilde{W}_{PSF}(\vec{k})$ and $\widetilde{W}_{\beta}(\vec{k})$
$[=\exp(-\beta^{2}\big{|}\vec{k}\big{|}^{2}/2)]$ are the Fourier transforms of
the PSF function and the Gaussian kernel respectively. $\beta$ is the scale
radius of the kernel, which is chosen to be slightly larger than that of the
original PSF, so that the reconvolution defined by $T(\vec{k})$ is
mathematically well defined. It can be shown that the FQ shear estimators
defined in eq.(2.1) can recover the underlying shear signals to the second
order in accuracy:

$\frac{<G_{1}>}{<N>}=g_{1}+\mathcal{O}(g_{1,2}^{3}),\;\frac{<G_{2}>}{<N>}=g_{2}+\mathcal{O}(g_{1,2}^{3})$
(2.3)

Although the form defined in eq.(2.3) is unbiased, it is not statistically
optimal, as the amplitudes of the $G_{i}$ and $N$ are unnormalized moments and
therefore strongly depend on the galaxy luminosity. Certain weightings can be
applied to reduce the scatter of the shear estimators, but the weighting
function should be carefully chosen to avoid additional biases [26].
Alternatively, in this work, we adopt the PDF-SYM method proposed by [23] to
measure the stacked shear signal and the shear-shear correlation. In this new
method, the shear signal is recovered by symmetrizing the probability
distribution function (PDF) of the FQ shear estimator $\hat{G_{i}}$ with an
assumed shear value $\hat{g_{i}}$ through
$\hat{G_{i}}=G_{i}-\hat{g_{i}}B_{i}$, in which $B_{1}=N+U$, $B_{2}=N-U$. Note
that the term $U$ here and $V$ defined in eq.(2.1) are the two components of a
spin-4 quantity, therefore both of them should be kept in the catalog.
Similarly, to measure the shear-shear correlation, the joint PDF of two shear
estimators $(\hat{G_{i}},\hat{G_{i}}^{\prime})$ is symmetrized using an
assumed correlation strength. With this method, we can bring the statistical
error down to the lower bound without introducing systematic errors. The
details of the PDF-SYM method can be found in [23]. Its performance has been
proved in a number of recent works [28, 29, 30, 31]. Its application on the
measurement of the shear-shear correlation is also discussed specifically in
our another recent work [32].

##  3 Tele-Correlation

For simplicity and clarity, we assume that the shear estimator takes a
conventional form $e_{i}$ with $i=1,2$, which can be regarded as the two
galaxy ellipticity components. Let us also denote the cosmic shear components
and the convergence as $\gamma_{i}$ and $\kappa$. The reduced shear $g_{i}$ is
therefore given by $\gamma_{i}/(1-\kappa)$. Under the weak lensing effect and
the field distortion effect, the shear estimator would take the following
form:

$e_{i}=e^{I}_{i}+(1+m_{i})(g_{i}+g^{f}_{i})+n_{i}$ (3.1)

where $e^{I}_{i}$ is the intrinsic ellipticity, and $g^{f}_{i}$ is the field
distortion signal. $m_{i}$ is the multiplicative bias associated with the
shear estimator, and $n_{i}$ is the noise, possibly containing an additive
bias. Tele-Correlation(TC) is defined as the correlation between the shear
estimators of two galaxies separated by a large distance (e.g., more than a
hundred degrees). In this case, the astrophysical correlation should vanish,
and TC only receives contribution from the field distortion:

$\displaystyle<e_{i}e_{i}^{\prime}>(\Delta\theta\rightarrow\infty)$
$\displaystyle=$
$\displaystyle(1+m_{i})(1+m_{i}^{\prime})<g_{i}^{f}g_{i}^{f\prime}>+<n_{i}n_{i}^{\prime}>$
$\displaystyle=$
$\displaystyle(1+\mathcal{M}_{i})^{2}<g_{i}^{f}g_{i}^{f\prime}>+\mathcal{C}_{i}$

We denote TC as $<e_{i}e_{i}^{\prime}>(\Delta\theta\rightarrow\infty)$ in the
rest of the paper. It is clear that by grouping the galaxy pairs according to
the products of their FD signals, TC can be directly compared with
$<g_{i}^{f}g_{i}^{f\prime}>$ to estimate $\mathcal{M}_{i}$ and
$\mathcal{C}_{i}$, and thereby the multiplicative and additive biases of the
shear estimators. This calibration can be done for galaxies of different
redshift bins. For the auto-correlation between the same population of the
galaxies, we should have $\mathcal{M}_{i}=m_{i}$.

Figure 1: TC of galaxy pairs with angular distance larger than 140 degrees. In
each panel, the black solid line is ’y=x’, and the red dashed line is the
best-fit line of data points. The $TC_{11}$ and $TC_{22}$ parts are shown in
the left and right panels respectively. The multiplicative and additive biases
for the left and right panels are $\mathcal{M}_{1}=-(1.99\pm 1.72)\times
10^{-2}$, $\mathcal{C}_{1}=(4.1\pm 2.6)\times 10^{-8}$,
$\mathcal{M}_{2}=(0.48\pm 0.78)\times 10^{-2}$, and $\mathcal{C}_{2}=(6.4\pm
1.53)\times 10^{-8}$ respectively.

Fig.1 shows the result of TC using the galaxy pairs separated by more than 140
degrees. We only choose the galaxies with redshifts larger than 0.1. The FD
product (x-axis in the figure) is evenly divided into 10 bins. The red dashed
lines are the best-fit of the data points, and the black solid lines refer to
’y=x’. According to the figure, the multiplicative biases are not significant:
$m_{1}=(-1.99\pm 1.717)\times 10^{-2}$; $m_{2}=(0.48\pm 0.779)\times 10^{-2}$,
and the additive biases $\mathcal{C}_{i}$ are negligible $(\sim 10^{-8})$.

It is interesting to note that one may expect to get smaller error bars in TC
if we enlarge the angular range of the galaxy pairs to get more samples. This
is indeed the case: we find that the average error-bar size of the data points
in fig.1 does decay when the lower limit of the separation angle is reduced.
However, the values and the uncertainties of the multiplicative biases do not
change much. In practice, when the lower limit of the galaxy angular
separation is set at 140 degree, the multiplicative biases converge well
enough, and the computational cost is already quite substantial777More than
70000 cores$\cdot$h for computing $1.7*10^{15}$ galaxies pairs.. As a
consistency check, we also set the angular range of the galaxy pairs to be
between 120 to 140 degrees, and the results are similar to those shown in
fig.1.

Figure 2: The TC of galaxies in different redshift bins using the photo-z from
[24]. $TC_{11}$ and $TC_{22}$ are shown in the upper and lower rows. In each
panel, the redshift range is shown in the upper-left corner. The red dashed
line is the best-fit for the data points, and the black solid line is the
’y=x’ line. The multiplicative and additive biases are shown in table.1.
Figure 3: The results of the original FD test (shear stacking) for different
redshift bins. In each panel, the redshift range is indicated in the upper-
left corner, and the red dashed line is the best-fit of the data points.
Figure 4: TC between different redshift bins as a consistency check. In each
panel, the label on the upper-left corner shows the two bin indices involved
in the measurement. The multiplicative bias $M$ is derived by fitting
$y=(1+M)^{2}x$ to the data points. The best-fit line is shown in red in each
panel. $M$ is expected to be consistent with the biases $m^{i}$ (”i” stands
for the index of the redshift bin) estimated from the auto-correlations, i.e.,
we expect:$(1+m^{i})(1+m^{j})=(1+M)^{2}$. For this purpose, we also show in
each panel the green line defined as $y=(1+m^{i})(1+m^{j})x$, using the bias
values from table.1. In this test, we only use the galaxy pairs with
separation larger than 170 degrees. Figure 5: The TC of galaxies in different
redshift bins. The photo-z catalog is from [25]. The definitions of the red
dashed lines and black solid lines are same as those in fig.2.

For the tomography study, we use TC to calibrate the shear biases in each
redshift bin. We divide the redshift range of $[0.2,1.0]$ evenly into 4 bins.
The TC results of these bins are shown in fig.2. The definitions of the axes
and the meanings of the black and red lines in the plots are the same as those
in fig.1. To our surprise, there is a significant multiplicative biases for
both $m_{1}$ and $m_{2}$ in a number of redshift bins. The values of the
biases are given in table.1. As a consistency check, we perform the original
FD tests (shear stacking) for the galaxies in these redshift bins. The results
are shown in fig.3 and table.2. The biases are consistent with those from the
TC in all the bins.

As a further consistency test, we also measure the TC between two different
redshift bins. The results are shown in fig.4. The label on the upper-left
corner of each panel show the indices of the two redshift bins. The bin index
of 1, 2, 3, 4 refers to the redshift range of $[0.2,0.4]$, $[0.4,0.6]$,
$[0.6,0.8]$, $[0.8,1]$ respectively. The definitions of the axes as well as
the black solid line and the red dashed line are the same as those in fig.2.
In this calculation, we only use the galaxy pairs separated by more than 170
degrees. Note that from the auto-correlation of each redshift bin we have
already obtained its multiplicative bias $m^{i}$ with i=1, 2, 3, 4 (note that
the upper index on $m$ refers to the redshift bin index). We expect that the
joint multiplicative bias $\mathcal{M}$ of the TC between two distinct bins is
given by $(1+\mathcal{M})^{2}=(1+m^{i})(1+m^{j})$. These predictions are shown
in fig.4 with the green lines given by $y=(1+\mathcal{M})^{2}x$. One can see
that the red and green lines in each plot agree well with each other, implying
that the multiplicative biases derived from the TC test are reliable.

It is still unclear to us why the redshift binning causes the multiplicative
biases, which can be as high as 10% or even more. There seem to be a strong
selection effect associated with the photo-z. If we switch to a different
photo-z catalog from [25], the biases persist, as shown in fig.5 and table.3.
Their amplitudes are similar to the original results in table.1.

We give more discussions in the §5 regarding the bias. Interestingly, in the
DES Y3 shear calibration simulation [15], they also find a redshift-dependent
multiplicative bias shown in different bins. We note that unlike the
conclusion of [15], the bias found in this work is NOT due to blending in our
test, because the blended sources experience the same FD effect. Given that
the FQ shear measurement does not need assumptions about the regularity of the
galaxy morphology, we do not expect shear bias to rise in the FD test.

In the next section, we implement these biases in our tomography study, and
check their influence on the cosmological constraints.

z | $10^{2}\times m_{1}$ | $10^{8}\times\mathcal{C}_{1}$ | $10^{2}\times m_{2}$ | $10^{8}\times\mathcal{C}_{2}$
---|---|---|---|---
[0.1-$\infty$] | -2.0$\pm$1.72 | 4.1$\pm$2.6 | 0.5$\pm$0.78 | 6.4$\pm$1.53
[0.2-0.4] | 17.5 $\pm$ 3.02 | -5.2$\pm$3.71 | 25.6 $\pm$ 3.26 | -3.6$\pm$5.26
[0.4-0.6] | -11.8 $\pm$ 5.34 | 7.0$\pm$6.69 | -10.9 $\pm$ 1.61 | 2.1$\pm$2.54
[0.6-0.8] | -3.2 $\pm$ 3.84 | -14.4$\pm$4.70 | -15.4 $\pm$ 2.70 | -3.3$\pm$4.04
[0.8-1.0] | -4.3 $\pm$ 8.72 | -8.5$\pm$11.30 | 12.9 $\pm$ 5.51 | 25.6$\pm$9.45
Table 1: The multiplicative and additive biases for different redshift bins from TC. The photo-z catalog is from [24]. z | $10^{2}\times m_{1}$ | $10^{4}\times c_{1}$ | $10^{2}\times m_{2}$ | $10^{4}\times c_{2}$
---|---|---|---|---
[0.1-$\infty$] | 3.2$\pm$4.36 | -0.52$\pm$0.35 | 0.9$\pm$4.18 | -0.44$\pm$0.35
[0.2-0.4] | 16.3 $\pm$ 3.64 | -0.8$\pm$0.29 | 24.9 $\pm$ 8.28 | -1.7$\pm$0.70
[0.4-0.6] | -6.6 $\pm$ 9.25 | -0.2$\pm$0.74 | -13.3 $\pm$ 3.57 | 0.6$\pm$0.30
[0.6-0.8] | -7.2 $\pm$ 6.71 | -0.7$\pm$0.54 | -14.7 $\pm$ 8.24 | 1.8$\pm$0.70
[0.8-1.0] | -9.7 $\pm$ 17.97 | 1.6$\pm$1.44 | 11.1 $\pm$ 26.44 | -4.1$\pm$2.24
Table 2: The multiplicative and additive biases for different redshift bins from the original FD test (shear stacking). The photo-z catalog is from [24]. z | $10^{2}\times m_{1}$ | $10^{8}\times\mathcal{C}_{1}$ | $10^{2}\times m_{2}$ | $10^{8}\times\mathcal{C}_{2}$
---|---|---|---|---
[0.2-0.4] | 13.6 $\pm$ 2.44 | -9.0$\pm$3.02 | 24.6 $\pm$ 2.90 | 0.1$\pm$4.87
[0.4-0.6] | -8.0 $\pm$ 6.11 | 1.7$\pm$7.27 | -8.5 $\pm$ 1.59 | 4.9$\pm$2.53
[0.6-0.8] | 1.3 $\pm$ 2.93 | 3.2$\pm$3.74 | -8.5 $\pm$ 2.16 | -0.8$\pm$3.49
[0.8-1.0] | -21.6 $\pm$ 7.13 | 16.2$\pm$10.05 | 1.3 $\pm$ 11.01 | 23.6$\pm$17.44
Table 3: The multiplicative and additive biases for the different redshift
bins from TC, using the photo-z catalog from [25].

##  4 Impact on Cosmological Constraints

###  4.1 Theory

The conventional shear-shear correlations are measured using the tangential
and cross shear components ($\gamma_{t}$ and $\gamma_{\times}$), which are
defined as $\gamma_{t}+i\gamma_{\times}=-(\gamma_{1}+i\gamma_{2})e^{-2\alpha
i},\alpha$ is the angle between the x-axis and the line connecting the galaxy
pair. Correspondingly, the correlations are defined as:

$\displaystyle\xi_{\pm}(z_{1},z_{2},\Delta\vec{\theta})$
$\displaystyle=<\gamma_{t}(z_{1},\vec{\theta})\gamma_{t}(z_{2},\vec{\theta}+\Delta{\vec{\theta}})>$
(4.1)
$\displaystyle\pm<\gamma_{\times}(z_{1},\vec{\theta})\gamma_{\times}(z_{2},\vec{\theta}+\Delta{\vec{\theta}})>$

In this paper, similar to the work of [32], we also consider another way:
correlations with only $\gamma_{1}$ or $\gamma_{2}$, i.e.,
$\xi_{ii}=<\gamma_{i}(z_{1},\vec{\theta})\gamma_{i}(z_{2},\vec{\theta}+\Delta{\vec{\theta}})>$,
which i =1 or 2. Theoretically, we expect $\xi_{11}=\xi_{22}=\xi_{+}/2$. We
use these different types of correlations to constrain cosmology.

In the theoretical model of $\xi_{+/-}$, one needs to consider the intrinsic
alignment. For our purpose, we simply adopt the recipe of [1], in which the
shear-shear correlation can be written as:

$\xi_{\pm}=\xi_{GG\pm}+\xi_{GI\pm}+\xi_{II\pm}$ (4.2)

where the subindex ”G” represents the cosmic shear, and ”I” stands for the
intrinsic galaxy shape. $\xi_{+/-}$ for galaxy pairs with the angular distance
$\theta$ are given by:

$\displaystyle\xi_{\pm}^{ij}(\theta)=\int_{0}^{\infty}\frac{d\ell}{2\pi}\ell
C^{ij}(\ell)J_{v}(\ell\theta),$ (4.3) $\displaystyle
C^{ij}(\ell)=\int_{0}^{\infty}dz\frac{W^{i}(z)W^{j}(z)}{\chi(z)^{2}}P_{\delta}(\frac{\ell}{\chi(z)},z)$

where $J_{\nu}$ is the first kind Bessel function, and $\nu$ is 0/4 for
$\xi_{+}$/$\xi_{-}$. $\chi(z)$ is the comoving radial distance at z.
$W^{i}(z)$ represents the kernel including the contributions from both lensing
and the intrinsic alignment, i.e.,
$W^{i}(z)=W_{\mathrm{G}}^{i}(z)+W_{I}^{i}(z)$. The kernel for shear is given
by:

$W_{\mathrm{G}}^{i}(z)=\frac{3}{2}\Omega_{m}\frac{H_{0}^{2}}{c^{2}}\frac{\chi(z)}{a(z)}\int_{z}^{\infty}dz^{\prime}n^{i}(z^{\prime})\frac{\chi(z^{\prime})-\chi(z)}{\chi(z^{\prime})},$
(4.4)

in which $n^{i}(z)$ is the normalized redshift distributions, $a(z)$ is the
scale factor, $H_{0}$ is the Hubble constant, and c is the speed of light. For
simplicity we assume a flat universe in this paper. For the IA kernel we adopt
the Nonlinear Alignment Model (NLA) [33]:

$W_{I}^{i}(z)=-\frac{A_{IA}C_{1}\rho_{c}\Omega_{m}}{D(z)}n^{i}(z)$ (4.5)

where $A_{IA}$ is a free parameter to show the amplitude of IA. $\rho_{c}$ is
the critical density, D(z) is the normalized growth factor. The normalization
constant $C_{1}$ is adjustable and can be set to $5\times
10^{-14}h^{-2}\mathrm{M}_{\odot}^{-1}\mathrm{Mpc}^{3}$ to align with the
observational findings[34]. In this case, the fiducial value of $A_{IA}$ is 1.

$P_{\delta}$ in eq.(4.3) is the nonlinear matter power spectrum, which could
be strongly affected by the baryonic effect. In our work, we use the baryonic
correction model (BCM, [35]) which parameterizes the influence of gas and
stars on the total matter density. There are two crucial parameters of the BCM
model: the mass fraction of ejected gas ($M_{c}$) and the ejection radius
(which depends on the parameter $\eta_{b}$). We choose their fiducial values
to be: $M_{c}=1.2M_{\odot}/h,\eta_{b}=0.5$. These values are consistent with
simulations and observations [35]. Overall, our theoretical predictions are
calculated by the Core Cosmology Library (CCL, [36]), in which the nonlinear
evolution is described by the halofit model [37, 38].

###  4.2 Redshift Distribution

In our tomography study, we use the photo-z catalog from [24]. For each
redshift bin selected according to the photometric redshifts, we hope to get
its true redshift distribution $n^{i}(z)$ used in eq.(4.3). It is related to
the photo-z distribution $f_{p}(z_{p})$ through the following equation:

$n^{i}(z_{s})=\int_{z^{i}_{min}}^{z^{i}_{max}}P(z_{s}|z_{p})f_{p}(z_{p})dz_{p}$
(4.6)

in which $z^{i}_{min}$ and $z^{i}_{max}$ are the lower and upper bounds of the
bin, $f_{p}(z_{p})$ is the overall photo-z distribution, and $P(z_{s}|z_{p})$
is the probability that a galaxy of $z_{p}$ has a true redshift of $z_{s}$.
Unfortunately, $P(z_{s}|z_{p})$ is not directly known to us. Its counterpart,
$P(z_{p}|z_{s})$, is usually more directly accessible by studying the
performance of photo-z reconstruction using simulations. They are related
through the Bayes theorem:

$P(z_{s}|z_{p})f_{p}(z_{p})=P(z_{p}|z_{s})f_{s}(z_{s}),$ (4.7)

where $f_{s}(z_{s})$ is the distribution of the true redshift, which is
neither known. If we can get an approximate form of $f_{s}(z_{s})$, the form
of $P(z_{s}|z_{p})$ can then be derived. For this purpose, we assume
$P(z_{p}|z_{s})$ is a Gaussian function with a pre-determined scatter
$\sigma_{z}$. Integrating both sides of eq.(4.7) over all possible values of
$z_{s}$, we get:

$f_{p}(z_{p})=\int f_{s}(z_{s})P(z_{p}|z_{s})dz_{s}$ (4.8)

In principle, $f_{s}(z_{s})$ can be derived by inverting the convolution in
eq.(4.8). In practice, we find that the solution from such an inversion is not
very stable. Assuming $\sigma_{z}$ is small, the integration in eq.(4.8) can
be well approximated as:

$f_{p}(z_{p})\approx f_{s}(z_{p})\int
P(z_{p}|z_{s})dz_{s}=f_{s}(z_{p})*g(z_{p})$ (4.9)

To get the form of $g(z)$, one can further convolve $f_{p}(z)$ with the same
Gaussian kernel $P(z_{p}|z_{s})$ to get $f_{c}(z)$, and therefore
$g(z_{p})\approx f_{c}(z)/f_{p}(z)$. Consequently, we get $f_{s}(z)\approx
f_{p}(z)/g(z)\approx f_{p}^{2}(z)/f_{c}(z)$. In the left panel of fig.6, we
show using simulations how well the $f_{s}(z)$ derived using the above method
can recover the true redshift distribution. The right panel of the same figure
shows the recovered redshift distribution for the four photo-z bins. We have
assumed that $P(z_{p}|z_{s})$ is a Gaussian kernel with
$\sigma_{z}=0.03*(1+z)$. This is a reasonable choice according to both [24]
and our another work (Li et al., in preparation). This assumption is used in
our tomography study in the next section. We have checked that our main
results do not change much if we assign $\sigma_{z}=0.05*(1+z)$ instead.

Figure 6: The left panel shows how well the true redshift distribution
$f_{s}(z)$ can be recovered with $f_{p}^{2}(z)/f_{c}(z)$ using the method
described in §4.2. The y-axis of its upper panel has an arbitrary unit, and
that of the lower panel shows the relative differences between $f_{p}(z)$ &
$f_{s}(z)$ (in blue), and $f_{p}^{2}(z)/f_{c}(z)$ & $f_{s}(z)$ (in red). The
right panel shows the redshift distributions of the four photo-z bins used for
our tomography study.

###  4.3 Tomography

Figure 7: The tomographic shear-shear correlation functions (2PCF). The
indices of the redshift bins are shown in the upper-left corner of each panel.
The blue and orange data points show the 2PCF without and with the corrections
from TC respectively. The blue points are shifted to their left by 10% for the
clearness of view. The black theoretical curves show the prdictions from [5]
with $A_{IA}=1$. The red curves are from the best-fit cosmological parameters
of the orange data points. The data in the shaded areas are excluded in our
cosmological constraints.

Fig.7 shows the measured two-point correlation functions
$\xi^{ij}_{\pm}(\theta)$ for all the pairs of redshift bins. The angular range
is from 1 arcmin to 300 arcmin. For the cosmological constraint, we do not use
the data points below 4 arcmin or above 100 arcmin. The blue points are the
results of $\xi_{\pm}$ without any calibration. The orange data points are the
results of $\xi_{\pm}$ after we incorporate the corrections from TC, i.e., we
multiply the shear estimator $G$ of each galaxy by $(1+m)^{-1}$ from the TC
according to its redshift bin. For clarity in the figure, we move the blue
points slightly to the left of their original places. The error bars and the
covariance matrix of the data points are estimated using the Jackknife method.
We use the K-Means clustering method [39] to divide the galaxies into 200
groups. The red dashed line is the best-fit cosmological prediction to the
orange data points (after calibration). The black dashed lines are calculated
from the cosmological parameters of PLANCK [5], with $A_{IA}=1$.

For the cosmological constraint, we use the standard Markov-Chain Monte Carlo
(MCMC) method with EMCEE[40] for cosmological parameter estimation. We treat
$A_{IA}$ of the intrinsic alignment as a free parameter. Our main results are
shown in fig.8. The blue contours and the red ones are for the best-fit
cosmological parameters from the blue and orange data points in fig.7
respectively. Without correcting for the redshift-dependent bias, we get:
$S_{8}=0.760^{+0.015}_{-0.017},\Omega_{m}=0.250^{+0.052}_{-0.037},A_{IA}=1.224^{+0.203}_{-0.203}$.
After we incorporate the corrections of the biases estimated from TC in each
redshift bin, we get:
$S_{8}=0.777^{+0.016}_{-0.019},\Omega_{m}=0.291^{+0.055}_{-0.048},A_{IA}=0.638^{+0.234}_{-0.264}$.
It is interesting and perhaps important to note that the $S_{8}$ value has
about 1$\sigma$ increase as a result of the calibration from TC.

Figure 8: Constraints on $S_{8}$, $\Omega_{m}$, and $A_{IA}$ using either
$\xi_{+}$ & $\xi_{-}$ or $\xi_{11}$ & $\xi_{22}$. The error contours in each
panel and each color correspond to the 68%, 95% and 99% confidence level
respectively. The blue and red contours show the constraints without and with
the corrections from TC respectively.

For a comparison, we show the results of $2*\theta\xi_{11}$ and
$2*\theta\xi_{22}$ in fig.9. The definitions of the blue and orange data
points, as well as the black and red lines are the same as those in fig.7. The
cosmological constraint using $\xi_{11}$ and $\xi_{22}$ together is shown in
the right panel of fig.8. Similarly, we find that the calibration from TC
leads to a $1\sigma$ increase in $S_{8}$. It is worth noting that the data
quality of $\xi_{11}$ and $\xi_{22}$ is somewhat worse than that of $\xi_{+}$
and $\xi_{-}$. The cosmological constraints from $\xi_{11}$ or $\xi_{22}$
separately are shown in fig.10, from which one can see that the quality of
$\xi_{11}$ is even worse than $\xi_{22}$. Similar phenomenon has been reported
in our another work [32].

In the DECaLS shear catalog, the first and second shear components are defined
in the local coordinates along the directions of ”RA” and ”Dec”, and the CCDs
are always lined up with the same direction in the survey, the image quality
therefore naturally inherits certain anisotropy from hardware imperfectness.
This problem is mitigated by rotating the shear estimators in the measurement
of $\xi_{+/-}$, yielding results of seemingly higher qualities. However, we
shall instead take this as a caution for the existence of unresolved
systematic issues in image processing.

Figure 9: Similar to fig.7, except that the cosmological constraints are from
$\xi_{11}$ & $\xi_{22}$. Figure 10: Similar to fig.8, but for the constraints
from $\xi_{11}$ (left) and $\xi_{22}$ (right) separately. Figure 11: Summary
of the $S_{8}$ constraints from various consistency tests, and from other
lensing surveys and Planck. Figure 12: Comparison of the TC results from
galaxies of different PSF sizes. In each panel, the red and blue points are
from galaxies with smaller (${\mathrm{FWHM_{PSF}}}<1.4"$) and larger
(${\mathrm{FWHM_{PSF}}}>1.4"$) PSF sizes respectively. The dashed lines are
the best-fit of the data points. The redshift range of the galaxies is
indicated in the upper-left corner of each panel. Figure 13: TC of galaxies in
the redshift range of [0.4, $\infty$]. The black solid lines and the dashed
lines have the same meanings as those in fig.1.

##  5 Conclusion and Discussions

Tele-correlation (TC), the correlation of the shear estimators of two galaxies
separated by a large distance ($\gtrsim$ 100 degree), can be used to calibrate
the multiplicative and additive biases in shear-shear correlations. In TC, the
correlation signal comes from the field distortion (FD), which can be retained
in the shear catalog if the shear estimators are associated with individual
exposures. We demonstrate this idea with the DECaLS shear catalog produced by
the Fourier_Quad pipeline. With all the distant galaxy pairs ( $>140$
degrees), we do not observe any bias in the correlation signal at the level of
$10^{-6}$. However, to our surprise, significant multiplicative biases can
arise if the TC test is performed in individual photo-z bins. The reason is
still not known to us.

Using the conventional tomographic shear-shear correlations $\xi_{+/-}$, we
try to place constraints on the cosmological parameters, and meanwhile study
the impact of the biases (from TC) on the results. We use the galaxies in the
redshift range of [0.2, 1.0], and the angular range from 4 to 100 arcmins for
calculating the correlations. We choose the angular range conservatively to
avoid the effects of some unknown physics at very small and large scales. In
our cosmological model, we adopt the NLA model for the intrinsic alignment,
and BCM for the impact of baryons on the density power spectrum. Without the
correction for the biases, we obtain $S_{8}=0.760^{+0.015}_{-0.017}$,
$\Omega_{m}=0.250^{+0.052}_{-0.037}$ and $A_{IA}=1.224\pm 0.203$. After
incorporating the multiplicative biases due to redshift binning, the
cosmological constraints become:
$S_{8}=0.777^{+0.016}_{-0.019},\Omega_{m}=0.291^{+0.055}_{-0.048},A_{IA}=0.638^{+0.234}_{-0.264}$.
There is about a $1\sigma$ increase in the best-fit value of $S_{8}$.

In fig.11, we summarize our constraints on $S_{8}$ from different types of
shear-shear correlations, including $\xi_{ii}$ ($=<\gamma_{i}\gamma_{i}>$, i=1
or 2), separately or jointly, with or without calibration from TC. In these
cases, we see a large variation of the constraints, particularly those with
$\xi_{11}$. This is an indication of the potential systematic biases in the
shear catalog, likely due to the image quality, as discussed in [32]. We also
present the $S_{8}$ constraints from the other lensing surveys such as the
fiducial $\Lambda$CDM-optimized analyses in DES Y3 [41], KiDS-1000 cosmology
[4], the analyses of $C_{l}$ and $\xi_{\pm}$ in HSC Y3 [3]. The Planck results
[5] with the baseline TT, TE, EE+LowE are also included. Our tomographic
lensing constraints on $S_{8}$ are mostly consistent with those of the other
lensing surveys, having about $2-3\sigma$ tension with the Planck results.

Unfortunately, we do not yet have a good understanding of the significant
multiplicative biases found in our TC test. It is likely a selection effect
because the galaxy redshift affects multiple image properties, including the
apparent magnitude, size, shape (after convolving with the PSF). The origin of
these biases remains elusive, and we will investigate the underlying
mechanisms in a future work. For now, to go a little further in this topic, we
divide the galaxy sample into two groups based on the PSF size. We measure the
tele-correlations for the first and second shear components ($TC_{11}$ and
$TC_{22}$), and show the results in fig.12. The black solid line in each panel
is the ’y=x’ line. The red lines are from the galaxies with smaller PSF size
(${\mathrm{FWHM_{PSF}}}<1.4"$), and the blue ones for larger PSF
(${\mathrm{FWHM_{PSF}}}>1.4"$). It is interesting to note that the
multiplicative biases seem to decrease with the PSF size. This fact agrees
with our intuition: shape measurement is more sensitive to the galaxy
size/redshift when the spatial resolution of the image is poorer, i.e., when
the PSF size is larger, thereby leading to a larger selection bias. This is a
good news for those surveys with small PSFs, such as HSC (Liu et al., in
preparation).

Some previous works (e.g., [29, 28, 42]) have used the FQ shear catalog of
DECaLS to study a number of statistics mainly in galaxy-galaxy lensing. We
realize that in those studies, since they only use the background galaxies
with redshifts larger than certain threshold (e.g., the lens redshift plus
some given value), the selection in redshift necessarily cause some bias that
needs to be corrected for. For example, we test the TC of the background
sample within the redshift range of $[0.4,\infty]$. The results are presented
in fig.13, which show a negative multiplicative bias of about $10\%$. Note
that this result can also be obtained with the original FD test. The actual
correction should be measured more carefully, as galaxies of different
redshifts are usually weighted by the critical surface density in galaxy-
galaxy lensing. The details will be discussed in our future work.

On the other hand, it is encouraging to note that the FD test offers us a
convenient onsite calibration of the shear biases, no matter whether the bias
is caused by PSF inaccuracy, instrumental effects, or selection effects. To
our knowledge, the only case in which the FD test may not work is about the
shear bias due to image blending, i.e., galaxy pairs that are close in angular
space, but far in redshift space. Since in this case, the two galaxies cannot
be treated as a single object in shear measurement, and the FD cannot tell
apart the two objects.

## Acknowledgments

This work is supported by the National Key Basic Research and Development
Program of China (2023YFA1607800, 2023YFA1607802), the NSFC grants (11621303,
11890691, 12073017), and the science research grants from China Manned Space
Project (No.CMS-CSST-2021-A01). The computations in this paper were run on the
Siyuan cluster supported by the Center for High-Performance Computing at
Shanghai Jiao Tong University.

## References

  * [1] C. Heymans, E. Grocutt, A. Heavens, M. Kilbinger, T.D. Kitching, F. Simpson et al., _CFHTLenS tomographic weak lensing cosmological parameter constraints: Mitigating the impact of intrinsic galaxy alignments_ , _MNRAS_ 432 (2013) 2433 [1303.1808].
  * [2] T.M.C. Abbott, M. Aguena, A. Alarcon, S. Allam, O. Alves, A. Amon et al., _Dark Energy Survey Year 3 results: Cosmological constraints from galaxy clustering and weak lensing_ , _Phys. Rev. D_ 105 (2022) 023520 [2105.13549].
  * [3] X. Li, T. Zhang, S. Sugiyama, R. Dalal, R. Terasawa, M.M. Rau et al., _Hyper Suprime-Cam Year 3 results: Cosmology from cosmic shear two-point correlation functions_ , _Phys. Rev. D_ 108 (2023) 123518 [2304.00702].
  * [4] M. Asgari, C.-A. Lin, B. Joachimi, B. Giblin, C. Heymans, H. Hildebrandt et al., _KiDS-1000 cosmology: Cosmic shear constraints and comparison between two point statistics_ , _A &A_ 645 (2021) A104 [2007.15633].
  * [5] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi et al., _Planck 2018 results. VI. Cosmological parameters_ , _A &A_ 641 (2020) A6 [1807.06209].
  * [6] R. Laureijs, J. Amiaux, S. Arduini, J.L. Auguères, J. Brinchmann, R. Cole et al., _Euclid Definition Study Report_ , _arXiv e-prints_ (2011) arXiv:1110.3193 [1110.3193].
  * [7] Y. Gong, X. Liu, Y. Cao, X. Chen, Z. Fan, R. Li et al., _Cosmology from the Chinese Space Station Optical Survey (CSS-OS)_ , _ApJ_ 883 (2019) 203 [1901.04634].
  * [8] R.T. Zellem, B. Nemati, V.P. Bailey, E.J. Cady, M.M. Colavita, G. Gonzalez et al., _Nancy Grace Roman Space Telescope Coronagraph Instrument Observation Calibration Plan_ , _arXiv e-prints_ (2022) arXiv:2202.05923 [2202.05923].
  * [9] M. Yamamoto, M.A. Troxel, M. Jarvis, R. Mandelbaum, C. Hirata, H. Long et al., _Weak Gravitational Lensing Shear Estimation with Metacalibration for the Roman High-Latitude Imaging Survey_ , _arXiv e-prints_ (2022) arXiv:2203.08845 [2203.08845].
  * [10] C. Heymans, L. Van Waerbeke, L. Miller, T. Erben, H. Hildebrandt, H. Hoekstra et al., _CFHTLenS: the Canada-France-Hawaii Telescope Lensing Survey_ , _MNRAS_ 427 (2012) 146 [1210.0032].
  * [11] L. Miller, C. Heymans, T.D. Kitching, L. van Waerbeke, T. Erben, H. Hildebrandt et al., _Bayesian galaxy shape measurement for weak lensing surveys - III. Application to the Canada-France-Hawaii Telescope Lensing Survey_ , _MNRAS_ 429 (2013) 2858 [1210.8201].
  * [12] X. Li, H. Miyatake, W. Luo, S. More, M. Oguri, T. Hamana et al., _The three-year shear catalog of the Subaru Hyper Suprime-Cam SSP Survey_ , _PASJ_ 74 (2022) 421 [2107.00136].
  * [13] A. Kannawadi, H. Hoekstra, L. Miller, M. Viola, I. Fenech Conti, R. Herbonnet et al., _Towards emulating cosmic shear data: revisiting the calibration of the shear measurements for the Kilo-Degree Survey_ , _A &A_ 624 (2019) A92 [1812.03983].
  * [14] M. Gatti, E. Sheldon, A. Amon, M. Becker, M. Troxel, A. Choi et al., _Dark energy survey year 3 results: weak lensing shape catalogue_ , _MNRAS_ 504 (2021) 4312 [2011.03408].
  * [15] N. MacCrann, M.R. Becker, J. McCullough, A. Amon, D. Gruen, M. Jarvis et al., _Dark Energy Survey Y3 results: blending shear and redshift biases in image simulations_ , _MNRAS_ 509 (2022) 3371 [2012.08567].
  * [16] J. Zhang, F. Dong, H. Li, X. Li, Y. Li, D. Liu et al., _Testing Shear Recovery with Field Distortion_ , _ApJ_ 875 (2019) 48 [1808.02593].
  * [17] J. Zhang, C. Liu, P.A. Vaquero, H. Li, H. Wang, Z. Shen et al., _Shear Measurement with Poorly Resolved Images_ , _AJ_ 164 (2022) 128 [2206.02434].
  * [18] H. Wang, J. Zhang, H. Li and Z. Shen, _Avoiding the geometric boundary effect in shear measurement_ , _The Astrophysical Journal_ 911 (2021) 10.
  * [19] J. Zhang, _Measuring the cosmic shear in Fourier space_ , _MNRAS_ 383 (2008) 113 [astro-ph/0612146].
  * [20] J. Zhang, _Measuring the reduced shear_ , _J. Cosmology Astropart. Phys._ 2011 (2011) 041 [1002.3614].
  * [21] J. Zhang and E. Komatsu, _Cosmic shears should not be measured in conventional ways_ , _MNRAS_ 414 (2011) 1047 [1002.3615].
  * [22] J. Zhang, W. Luo and S. Foucaud, _Accurate shear measurement with faint sources_ , _J. Cosmology Astropart. Phys._ 2015 (2015) 024 [1312.5514].
  * [23] J. Zhang, P. Zhang and W. Luo, _Approaching the Cramér-Rao Bound in Weak Lensing with PDF Symmetrization_ , _ApJ_ 834 (2017) 8 [1606.09397].
  * [24] R. Zhou, J.A. Newman, Y.-Y. Mao, A. Meisner, J. Moustakas, A.D. Myers et al., _The clustering of DESI-like luminous red galaxies using photometric redshifts_ , _MNRAS_ 501 (2021) 3309 [2001.06018].
  * [25] H. Zou, J. Gao, X. Zhou and X. Kong, _Photometric Redshifts and Stellar Masses for Galaxies from the DESI Legacy Imaging Surveys_ , _ApJS_ 242 (2019) 8.
  * [26] H. Li, J. Zhang, D. Liu, W. Luo, J. Zhang, F. Dong et al., _Toward a Bias-free Selection Criterion in Shear Measurement_ , _ApJ_ 908 (2021) 93 [2006.02095].
  * [27] T. Lu, J. Zhang, F. Dong, Y. Li, D. Liu, L. Fu et al., _Testing PSF Interpolation in Weak Lensing with Real Data_ , _AJ_ 153 (2017) 197 [1610.09828].
  * [28] M. Fong, J. Han, J. Zhang, X. Yang, H. Gao, J. Wang et al., _First measurement of the characteristic depletion radius of dark matter haloes from weak lensing_ , _MNRAS_ 513 (2022) 4754 [2205.01816].
  * [29] J. Wang, X. Yang, J. Zhang, H. Li, M. Fong, H. Xu et al., _Halo Properties and Mass Functions of Groups/Clusters from the DESI Legacy Imaging Surveys DR9_ , _ApJ_ 936 (2022) 161 [2207.12771].
  * [30] Z. Zhang, H. Wang, W. Luo, J. Zhang, H.J. Mo, Y. Jing et al., _Massive star-forming galaxies have converted most of their halo gas into stars_ , _A &A_ 663 (2022) A85 [2112.04777].
  * [31] H. Wang, J. Zhang, H. Li and C. Liu, _Toward an Optimal Reconstruction of the Shear Field with PDF-folding_ , _ApJ_ 954 (2023) 193 [2303.08313].
  * [32] Z. Liu, J. Zhang, H. Li, Z. Shen and C. Liu, _Quasi-2D weak lensing cosmological constraints using the PDF-SYM method_ , _Science China Physics, Mechanics, and Astronomy_ 67 (2024) 270413 [2310.11209].
  * [33] S. Bridle and L. King, _Dark energy constraints from cosmic shear power spectra: impact of intrinsic alignments on photometric redshift requirements_ , _New Journal of Physics_ 9 (2007) 444 [0705.0166].
  * [34] M.L. Brown, A.N. Taylor, N.C. Hambly and S. Dye, _Measurement of intrinsic alignments in galaxy ellipticities_ , _MNRAS_ 333 (2002) 501 [astro-ph/0009499].
  * [35] A. Schneider and R. Teyssier, _A new method to quantify the effects of baryons on the matter power spectrum_ , _J. Cosmology Astropart. Phys._ 2015 (2015) 049 [1510.06034].
  * [36] N.E. Chisari, D. Alonso, E. Krause, C.D. Leonard, P. Bull, J. Neveu et al., _Core Cosmology Library: Precision Cosmological Predictions for LSST_ , _ApJS_ 242 (2019) 2 [1812.05995].
  * [37] R.E. Smith, J.A. Peacock, A. Jenkins, S.D.M. White, C.S. Frenk, F.R. Pearce et al., _Stable clustering, the halo model and non-linear cosmological power spectra_ , _MNRAS_ 341 (2003) 1311 [astro-ph/0207664].
  * [38] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya and M. Oguri, _Revising the Halofit Model for the Nonlinear Matter Power Spectrum_ , _ApJ_ 761 (2012) 152 [1208.2701].
  * [39] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel et al., _Scikit-learn: Machine Learning in Python_ , _Journal of Machine Learning Research_ 12 (2011) 2825 [1201.0490].
  * [40] D. Foreman-Mackey, D.W. Hogg, D. Lang and J. Goodman, _emcee: The MCMC Hammer_ , _PASP_ 125 (2013) 306 [1202.3665].
  * [41] A. Amon, D. Gruen, M.A. Troxel, N. MacCrann, S. Dodelson, A. Choi et al., _Dark Energy Survey Year 3 results: Cosmology from cosmic shear and robustness to data calibration_ , _Phys. Rev. D_ 105 (2022) 023514 [2105.13543].
  * [42] H. Xu, H. Li, J. Zhang, X. Yang, P. Zhang, M. He et al., _DESI Legacy Imaging Surveys Data Release 9: Cosmological constraints from galaxy clustering and weak lensing using the minimal bias model_ , _Science China Physics, Mechanics, and Astronomy_ 66 (2023) 129811 [2310.03066].


#  Stagewise Newton Method for Dynamic Game Control
with Imperfect State Observation

Armand Jordana1, Bilal Hammoud1, Justin Carpentier2 and Ludovic Righetti1,3
1Tandon School of Engineering, New York University, Brooklyn, NY
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Département
d’informatique de l’ENS, École normale supérieure, CNRS, PSL Research
University, Paris, France<EMAIL_ADDRESS>Institute for
Intelligent Systems, Tübingen, GermanyThis work was in part supported by the
European Union Horizon 2020 research and innovation program (grant agreement
780684), the National Science Foundation grants 1825993, 1932187, 1925079,
2026479, and the French government under management of Agence Nationale de la
Recherche as part of the ”Investissements d’avenir” program, reference
ANR-19-P3IA-0001 (PRAIRIE 3IA Institute), Louis Vuitton ENS Chair on
Artificial Intelligence.

###### Abstract

In this letter, we study dynamic game optimal control with imperfect state
observations and introduce an iterative method to find a local Nash
equilibrium. The algorithm consists of an iterative procedure combining a
backward recursion similar to minimax differential dynamic programming and a
forward recursion resembling a risk-sensitive Kalman smoother. A coupling
equation renders the resulting control dependent on the estimation. In the
end, the algorithm is equivalent to a Newton step but has linear complexity in
the time horizon length. Furthermore, a merit function and a line search
procedure are introduced to guarantee convergence of the iterative scheme. The
resulting controller reasons about uncertainty by planning for the worst case
disturbances. Lastly, the low computational cost of the proposed algorithm
makes it a promising method to do output-feedback model predictive control on
complex systems at high frequency. Numerical simulations on realistic robotic
problems illustrate the risk-sensitive behavior of the resulting controller.

##  I Introduction

Doing output-feedback Model Predictive Control (MPC) while being robust to the
estimator uncertainty is a notoriously difficult problem [1]. In optimal
control, when the state measurement is partial and corrupted by noise, a
common practice is to treat the estimation problem independently from the
control and rely on the most likely outcome. However, in some scenarios, one
might want to design a controller robust to the worst case estimation error.
Such controllers can be obtained through dynamic game control [2, 3]. As Mayne
[1] advocates, one meaningful way to properly take into account the
measurement uncertainty is to use a minimax formulation linking control and
estimation.

This problem has been extensively studied in [4, 5], which shows that a
specific dynamic game formulation leads to MPC approaches with bounded state
trajectories and provides an explicit characterization of these bounds.
However, the minimax problem was solved with an interior point method without
taking into account the specific structure of the problem and the sparsity
induced by time. In this work, we derive an explicit iterative solution that
fully exploits sparsity, resulting in an algorithm that linearly scales with
the time horizon length and which can be easily warm-started for use in MPC
schemes [6].

In the linear dynamics and quadratic cost case, Jacobson [7] showed that
dynamic game control is equivalent to risk sensitive control and derived a
closed form solution. Later, Whittle [8] extended the results to the linear
quadratic case with imperfect state observations. In [9], a first iterative
version of Whittle’s solution was introduced to tackle the nonlinear risk
sensitive problem with imperfect observations. However, the stochastic nature
of the problem hindered the development of theoretical guarantees. Although,
dynamic game and risk sensitive control are equivalent in the linear quadratic
case [7], this is no longer the case in the nonlinear setting [10].
Nonetheless, dynamic game control is tightly connected to robust and risk
sensitive control. In [2, 10], James and Campi showed that dynamic game
control can be interpreted as the limit case of risk sensitive control when
noise tends to zero. Recently, Başar [11] presented a detailed overview of the
connections between both problems in continuous time. Additionally, Başar and
Bernhard [3, 12] established the connections between dynamic game control and
$\mbox{H}^{\infty}$-Optimal Control both in the perfect and imperfect state
information case.

For nonlinear systems, estimating a state trajectory corresponding to some
given measurements is usually intractable analytically. A common approach is
to model the noise as Gaussian and to maximize the Maximum A Posteriori (MAP)
[13]. If the dynamics are affine, then the problem can be solved analytically
with the so-called Rauch–Tung–Striebel (RTS) smoother [14]. The RTS smoother
is made of a forward recursion which resembles the Kalman Filter (KF) and a
backward recursion. In the nonlinear case, iterative schemes are usually used.
A popular choice is the iterative Kalman smoother which is equivalent to a
Gauss-Newton method on the MAP [15]. The estimation part of our proposed
solution resembles a risk sensitive version of this smoother.

In optimal control or dynamic game control with perfect state information,
various numerical optimization algorithms have been developed to iteratively
find solutions in the nonlinear case. In optimal control, the most analogous
to our work are Differential Dynamic Programming (DDP) [16] and the stagewise
implementation of the Newton’s method [17]. The stagewise Newton method is an
exact implementation of Newton’s method that exploits the specific structure
of the Hessian matrix in order to scale linearly with the time horizon. DDP is
an iterative algorithm that takes an update step on the control input by
applying dynamic programming on a quadratic approximation of the value
function. In [16], Murray showed that DDP is very similar to a Newton step and
inherits its convergence properties. For dynamic game control with perfect
state information, the seminal work from [18, 19] introduced minimax DDP
showing that DDP could be extended to zero-sum two players games. Recently,
[20] further extended the concepts of stagewise Newton method and DDP to
nonzero-sum games with an arbitrary number of players in the full information
case.

Despite having been widely studied theoretically, to the best of our
knowledge, dynamic game control with imperfect state observations has not been
approached from a numerical optimization point of view. In this work, we
consider the general problem of dynamic game control with imperfect state
observation and present a numerically efficient provably convergent algorithm
to solve it. The solution presented generalizes commonly used techniques such
as the extended Kalman smoother and Differential Dynamic Programming. The
proposed solution fully exploits the sparsity of the problem and scales
linearly with time, making it a promising method for model predictive control.

Contributions. In this paper, we derive a stagewise Newton method for dynamic
game control with imperfect state observations. The proposed approach scales
linearly in the time horizon. Furthermore, a merit function and a line search
procedure are introduced to guarantee convergence. To our knowledge, this is
the first iterative procedure with convergence guarantee for the dynamic game
control problem with imperfect state observation. The proposed solver couples
estimation and control by merging an iterative optimal control algorithm
similar to minimax DDP and an iterative risk sensitive Kalman smoother. We
illustrate the behavior of the resulting controller on two realistic robotic
problems.

Notations. The derivative of a function $f$ with respect to a vector $v$ is
denoted by $f^{v}$, similarly for second order derivatives with respect to
vectors $u,v$ is denoted as $f^{uv}$. If $(v_{i})_{i\in\mathbb{N}}$ is a
sequence of vectors, then $v_{k:t}$ is a vector concatenating $v_{k}\dots
v_{t}$. $\mathbf{1}_{x\in A}$ is the indicator function which equals $1$ if
$x\in A$ and $0$ otherwise. $I_{n}$ denotes the identity matrix of size $n$ by
$n$.

##  II Problem Statement

Similar to [4, 5], this work studies a special class of nonlinear dynamic
games with imperfect state observation [2]. Given a history of measurements
$y_{1:t}$, a history of control inputs, $u_{0:t-1}$ and a prior on the initial
state $\hat{x}_{0}$, we aim to find a control sequence $u_{t:T-1}$ that
minimizes a given cost $\ell$ while an opposing player aims to find the
disturbances $(w_{0:T},\gamma_{1:t})$ that maximize this cost $\ell$ minus a
weighed norm of the disturbances. Such a problem is formally written as:

$\displaystyle\min_{u_{t:T-1}}\max_{w_{0:T}}\max_{\gamma_{1:t}}\sum_{j=0}^{T-1}\ell_{j}(x_{j},u_{j})+\ell_{T}(x_{T})$
(1)
$\displaystyle-\frac{1}{2\mu}\left(\omega_{0}^{T}P^{-1}\omega_{0}+\sum_{j=1}^{t}\gamma_{j}^{T}R_{j}^{-1}\gamma_{j}+\sum_{j=1}^{T}w_{j}^{T}Q_{j}^{-1}w_{j}\right)$
$\displaystyle\mbox{subject to}\quad x_{0}$ $\displaystyle=\hat{x}_{0}+w_{0},$
(2a) $\displaystyle x_{j+1}$ $\displaystyle=f_{j}(x_{j},u_{j})+w_{j+1},$
$\displaystyle 0\leq j<T,$ (2b) $\displaystyle y_{j}$
$\displaystyle=h_{j}(x_{j})+\gamma_{j},$ $\displaystyle 1\leq j\leq t.$ (2c)

where $\mu>0$. $x_{j}$ is the state, $\omega_{j}$ the process disturbance,
$\gamma_{j}$ the measurement disturbance, $T$ the time horizon, $t$ the
current time. The transition model $f_{j}$, the measurement model $h_{j}$ and
the cost $\ell_{j}$ are assumed to be $\mathcal{C}^{2}$. The measurement
uncertainty $R_{j}$, the process uncertainty $Q_{j}$ and the initial state
uncertainty $P$ are positive-definite matrices.

Interestingly, this problem encompass various formulations of control and
estimation. If $t=0$ and if $w_{0}$ is fixed to zero, we recover dynamic game
control with perfect state information. Additionally, if $t=0$, in the limit
where $\mu$ tends to zero, we find the generic optimal control formulation
[10]. And lastly, if $t=T$ and if we consider all the cost $\ell_{j}$ to be
null, then, (1) is equivalent to maximizing the MAP.

##  III Method

###  III-A First order conditions for a Nash equilibrium

The main challenge in the problem formulation (1) is the equality constraint
maintaining the dynamic feasibility. A popular approach [19] is to derive a
DDP-like algorithm by sequentially taking quadratic approximations of the
value function recursion. However, a stagewise Newton step can readily be
derived. One of the key features of the dynamic game (1) is that by changing
the decision variable of the opposing player, one can transform the problem
into an unconstrained one [2]. Indeed, we can use the equality constraints of
Eq. (2) to replace disturbance maximization into a maximization over the state
sequence. Then the problem loses its equality constraints and can be
formulated as the search of the saddle point of

$\displaystyle J($ $\displaystyle
x_{0:T},u_{t:T-1})=\sum_{j=0}^{T-1}\ell_{j}(x_{j},u_{j})+\ell_{T}(x_{T})$ (3)
$\displaystyle-\frac{1}{2\mu}(x_{0}-\hat{x}_{0})^{T}P^{-1}(x_{0}-\hat{x}_{0})$
$\displaystyle-\frac{1}{2\mu}\sum_{j=1}^{t}(y_{j}-h_{j}(x_{j}))^{T}R_{j}^{-1}(y_{j}-h_{j}(x_{j}))$
$\displaystyle-\frac{1}{2\mu}\sum_{j=0}^{T-1}(x_{j+1}-f_{j}(x_{j},u_{j}))^{T}Q_{j+1}^{-1}(x_{j+1}-f_{j}(x_{j},u_{j})).$

However, without convexity and concavity assumptions, we cannot aim at finding
global solutions of the minimax problem. Hence, we restrict our attention to
local Nash equilibrium, namely a point $(x_{0:T}^{\star},u_{t:T-1}^{\star})$
such that there exits $\delta>0$ such that for any $(x_{0:T},u_{t:T-1})$
satisfying $||x_{0:T}-x^{\star}_{0:T}||<\delta$ and
$||u_{t:T-1}-u^{\star}_{t:T-1}||<\delta$, we have

$\displaystyle J(x_{0:T},u^{\star}_{t:T-1})\leq
J(x^{\star}_{0:T},u^{\star}_{t:T-1})\leq J(x^{\star}_{0:T},u_{t:T-1}).$ (4)

A standard approach to this problem is to search for a stationary point [21]:

$\displaystyle\begin{pmatrix}\frac{\partial J}{\partial
x_{0:T}}(x^{\star}_{0:T},u^{\star}_{t:T-1})\\\\[2.84526pt] \frac{\partial
J}{\partial u_{t:T-1}}(x^{\star}_{0:T},u^{\star}_{t:T-1})\end{pmatrix}=0$ (5)

Interestingly, the change of decision variable in Eq. 3 turned the problem
into an unconstrained one but made no assumption on the structure of the cost.
The only required assumption is that each disturbance is in a one to one map
with the state at each time step.

###  III-B About the Linear Quadratic case

In the linear quadratic case, Whittle has shown that under some conditions,
the saddle point of (4) is global and can be computed analytically. More
precisely, the order of the minimization and maximization in (1) can be
interchanged, namely the lower value and upper value of the game are equal.
Despite this result, one of the difficulties of the problem (1) is that it
links estimation and control. One the major contribution of Whittle is the
introduction of the notion of past stress and future stress showing that the
KF and LQR principles can still be applied. In other words, the problem can be
solved by performing a backward recursion on the controls and future states
and a forward recursion on the past states. The future stress recursion can be
interpreted as a value function recursion similar to LQR, while the past
stress can be interpreted as a rollout of KF. Here, we use these two
principles to efficiently solve iteratively the nonlinear case.

###  III-C A stagewise Newton’s method

In this section, a stagewise formulation of the Newton method to find a
stationary point of $J$ is introduced. While a naive implementation of the
Newton method would yield a complexity of $O(T^{3})$, it is shown that the
special structure of the Hessian induced by time can be exploited in order to
obtain a linear complexity in time $O(T)$. In the perfect state observation
case, [17] and [20] derived a stagewise Newton method with a backward
recursion on the controls. However, with imperfect state observation, it is no
longer clear how to do this with only one recursion. Instead, we show that the
principles introduced by Whittle can be applied.

To ensure that the proposed method is well defined and to guarantee
convergence, we assume that the cost satisfies smoothness and non-degeneracy
conditions required for the convergence of Newton’s method [22]. As the cost
(3) is unconstrained, the gradients and Hessian of the cost can readily be
computed. At iteration $i$, given a guess
$x_{0}^{i},x_{1}^{i}\dots,x_{T}^{i},u_{t}^{i}\dots,u_{T-1}^{i}$, also referred
as the nominal trajectory, the Newton step, denoted by $p$, satisfies

$\displaystyle\begin{pmatrix}\dfrac{\partial^{2}J}{\partial x_{0:T}\partial
x_{0:T}}&\dfrac{\partial^{2}J}{\partial x_{0:T}\partial
u_{t:T-1}}\\\\[8.53581pt] \dfrac{\partial^{2}J}{\partial u_{t:T-1}\partial
x_{0:T}}&\dfrac{\partial^{2}J}{\partial u_{t:T-1}\partial
u_{t:T-1}}\end{pmatrix}p=\begin{pmatrix}\dfrac{\partial J}{\partial
x_{0:T}}\\\\[8.53581pt] \dfrac{\partial J}{\partial u_{t:T-1}}\end{pmatrix}$
(6)

Here, $p=\begin{pmatrix}p_{x_{0:T}}^{T}&p_{u_{t:T-1}}^{T}\end{pmatrix}^{T}$
where ${p_{x_{0:T}}\in\mathbb{R}^{(T+1)n_{x}}}$ is a stack of vectors
$p_{x_{k}}\in\mathbb{R}^{n_{x}}$ with $n_{x}$ being the dimension of the state
space. Similarly, ${p_{u_{t:T-1}}\in\mathbb{R}^{(T-t)n_{u}}}$ is a stack of
vectors $p_{u_{k}}\in\mathbb{R}^{n_{u}}$ with $n_{u}$ being the dimension of
the control space. To simplify the notations, we define an augmented Hessian
of the cost that contains the second order derivatives of the dynamics for all
$k<T$.

$\displaystyle\bar{\ell}^{xx}_{k}$
$\displaystyle=\ell^{xx}_{k}+\mu^{-1}{w_{k+1}^{i}}^{T}Q_{k+1}^{-1}f^{xx}_{k}+\mu^{-1}\mathbf{1}_{1\leq
k\leq t}{\gamma_{k}^{i}}^{T}R_{k}^{-1}h^{xx}_{k},$
$\displaystyle\bar{\ell}^{xu}_{k}$
$\displaystyle=\bar{\ell}_{k}^{{ux}^{T}}=\ell^{xu}_{k}+\mu^{-1}{w_{k+1}^{i}}^{T}Q_{k+1}^{-1}f^{xu}_{k},$
$\displaystyle\bar{\ell}^{uu}_{k}$
$\displaystyle=\ell^{uu}_{k}+\mu^{-1}{w^{i}_{k+1}}^{T}Q_{k+1}^{-1}f^{uu}_{k},$
(7)

where the derivatives are evaluated at the current guess and where
${w_{k+1}^{i}\coloneqq x_{k+1}^{i}-f_{k}(x_{k}^{i},u_{k}^{i})}$ and
${\gamma_{k}^{i}\coloneqq y_{k}-h_{k}(x_{k}^{i})}$. Here, the second order
derivatives of the dynamics are tensors. The exact definition of the tensor
indexing and the tensor product is provided in the supplementary material
[23].

The next three propositions are analogous to the principles introduced by
Whittle: the past stress recursion, the future stress recursion and the
coupling of the past and future stress recursions. The first proposition,
analogous to the future stress recursion, expresses every future state and
control update steps as a function of $p_{x_{t}}$.

######  Proposition 1 (Future stress).

In Equation (6), the last ${(T-t)(n_{x}+n_{u})}$ rows are equivalent to:

$\displaystyle\forall k\geq t,\quad p_{u_{k}}=G_{k}p_{x_{k}}+g_{k}$ (8)
$\displaystyle p_{x_{k+1}}=\left(I-\mu Q_{k+1}V_{k+1}\right)^{-1}($
$\displaystyle f^{x}_{k}p_{x_{k}}+f^{u}_{k}p_{u_{k}}$ $\displaystyle+\mu
Q_{k+1}v_{k+1}-w^{i}_{k+1})$

where $V_{k}$ and $v_{k}$ are solutions of the backward recursion:

$\displaystyle\Gamma_{k+1}$ $\displaystyle=I-\mu V_{k+1}Q_{k+1}$ (9)
$\displaystyle Q_{uu}$
$\displaystyle=\bar{\ell}^{uu}_{k}+{f^{u}_{k}}^{T}\Gamma_{k+1}^{-1}V_{k+1}f^{u}_{k}$
$\displaystyle Q_{ux}$
$\displaystyle=\bar{\ell}^{ux}_{k}+{f^{u}_{k}}^{T}\Gamma_{k+1}^{-1}V_{k+1}f^{x}_{k}$
$\displaystyle Q_{u}$
$\displaystyle=\ell^{u}_{k}+{f^{u}_{k}}^{T}\Gamma_{k+1}^{-1}\left(v_{k+1}-V_{k+1}w^{i}_{k+1}\right)$
$\displaystyle G_{k}$ $\displaystyle=-Q_{uu}^{-1}Q_{ux}$ $\displaystyle g_{k}$
$\displaystyle=-Q_{uu}^{-1}Q_{u}$ $\displaystyle V_{k}$
$\displaystyle=\bar{\ell}^{xx}_{k}+{f^{x}_{k}}^{T}\Gamma_{k+1}^{-1}V_{k+1}f^{x}_{k}+Q_{ux}^{T}G_{k}$
$\displaystyle v_{k}$
$\displaystyle=\ell^{x}_{k}+{f^{x}_{k}}^{T}\Gamma_{k+1}^{-1}\left(v_{k+1}-V_{k+1}w^{i}_{k+1}\right)+Q_{ux}^{T}g_{k}$

with the terminal condition

$\displaystyle V_{T}$ $\displaystyle=\ell^{xx}_{T},$ $\displaystyle
v_{T}=\ell^{x}_{T}.$ (10)

In those equations, $\mu$ intervenes only in $\Gamma_{k}$ and the augmented
terms of the cost. Interestingly, $\Gamma^{-1}_{k}$ shifts, at each time step,
the value function terms $V_{k}$ and $v_{k}$. Then, the second proposition,
analogous to the past stress recursion, expresses every past state update
steps as a function of $p_{x_{t}}$.

######  Proposition 2 (Past stress).

In Equation (6), if $t\geq 1$, the first ${(t-1)n_{x}}$ rows are equivalent
to: $\forall k=0,\dots,t-1,$

$\displaystyle p_{x_{k}}$
$\displaystyle=E_{k+1}^{-1}\left({f^{x}_{k}}^{T}Q_{k+1}^{-1}(w_{k+1}^{i}+p_{x_{k+1}})+P_{k}^{-1}\hat{\mu}_{k}+\mu
l^{x}_{k}\right)$ (11)

where $P_{k}$ and $\hat{\mu}_{k}$ are solution of the forward recursion:

$\displaystyle E_{k+1}$
$\displaystyle=P_{k}^{-1}+{f^{x}_{k}}^{T}Q_{k+1}^{-1}f^{x}_{k}-\mu\bar{l}^{xx}_{k}$
$\displaystyle\bar{P}_{k+1}$
$\displaystyle=Q_{k+1}+f^{x}_{k}(P_{k}^{-1}-\mu\bar{\ell}^{xx}_{k})^{-1}{f^{x}_{k}}^{T}$
$\displaystyle K_{k+1}$
$\displaystyle=\bar{P}_{k+1}h_{k+1}^{x^{T}}(R_{k+1}+h^{x}_{k+1}\bar{P}_{k+1}h_{k+1}^{x^{T}})^{-1}$
$\displaystyle P_{k+1}$ $\displaystyle=(I-K_{k+1}h^{x}_{k+1})\bar{P}_{k+1}$
$\displaystyle\hat{\mu}_{k+1}$
$\displaystyle=(I-K_{k+1}h^{x}_{k+1})(f^{x}_{k}\hat{\mu}_{k}-w_{k+1}^{i})+K_{k+1}\gamma_{k+1}^{i}$
$\displaystyle+\mu
P_{k+1}Q_{k+1}^{-1}f^{x}_{k}E_{k+1}^{-1}(\bar{\ell}^{xx}_{k}\hat{\mu}_{k}+\ell^{x}_{k})$
(12)

with the initialization

$\displaystyle P_{0}$ $\displaystyle=P,$
$\displaystyle\hat{\mu}_{0}=\hat{x}_{0}-x_{0}^{i}.$ (13)

Interestingly, if all the cost terms $\ell_{j}$ are zero and if $t=T$, the
method is equivalent to a Newton method on the MAP. Furthermore, if the second
order derivatives of the measurement function are omitted, then, the algorithm
is equivalent to the iterative Kalman Smoother. Finally, the third proposition
shows how both past stress and future stress recursions can be coupled to find
the update step $p_{x_{t}}$.

######  Proposition 3 (Coupling).

In Equation (6), the remaining rows (from $tn_{x}+1$ to $(t+1)n_{x}$) are
equivalent to

$\displaystyle p_{x_{t}}=\left(P_{t}^{-1}-\mu
V_{t}\right)^{-1}\left(P_{t}^{-1}\hat{\mu}_{t}+\mu v_{t}\right).$ (14)

In the limit case when $\mu$ tends to zero, the estimation and control are
decoupled and we recover the usual certainty equivalence principle:
$p_{x_{t}}=\hat{\mu}_{t}$. The algorithm is then equivalent to an iterative
estimator and an iterative controller running independently. More precisely,
at each iteration the controller uses the current estimate of the smoother.

###### Proof.

The proof follows from the analytical derivations of the gradient and the
Hessian of (3), a forward induction from $0$ to $t$ and a backward induction
from $T$ to $t$. The complete proof is provided in the supplementary material
[23]. ∎

Input: $x_{0}^{i},x_{1}^{i}\dots,x_{T}^{i},u_{t}^{i}\dots,u_{T-1}^{i}$

// Estimation forward pass

1 $P_{0}\leftarrow P$, $\hat{\mu}_{0}\leftarrow\hat{x}_{0}-x_{0}^{i}$

2 for  _$k=0,...t-1$_ do

3 $P_{k+1},\hat{\mu}_{k+1}\leftarrow$ Eq. (12)

4

// Control backward pass

5 $V_{T}\leftarrow\ell^{xx}_{T}$, $v_{T}\leftarrow\ell^{x}_{T}$

6 for  _$k=T-1,...t$_ do

7 $V_{k},v_{k}\leftarrow$ Eq. (9)

8

// Estimation and control coupling

9 $p_{x_{t}}\leftarrow\left(P_{t}^{-1}-\mu
V_{t}\right)^{-1}\left(P_{t}^{-1}\hat{\mu}_{t}+\mu v_{t}\right)$

// Estimation backward pass

10 for  _$k=t-1,...0$_ do

11 $p_{x_{k}}\leftarrow$ Eq. (11)

12

// Control forward pass

13 for  _$k=t,\dots T-1$_ do

14 $p_{u_{k}},p_{x_{k+1}}\leftarrow$ Eq. (8)

15

Output: $p_{x_{0}},p_{x_{1}}\dots,p_{x_{T}},p_{u_{t}}\dots,p_{u_{T-1}}$

Algorithm 1 Stagewise Newton step

In the end, the update step, $p$, can be computed with a forward recursion on
the past indexes, a backward recursion on the future indexes, a coupling
equation, a backward recursion on the past indexes and a forward recursion on
the future indexes. Algorithm 1 summarizes those steps. Clearly, the
complexity is linear in time. Instead of inverting a matrix of size
${((T+1)n_{x}+(T-t)n_{u})}$, Algorithm 1 only operates with matrices of size
$n_{x}$ or $n_{u}$ and the number of operation is proportional to $T$.

###  III-D Line-search and convergence

In the linear quadratic case, Algorithm 1 is equivalent to Whittle’s
derivations and only one iteration is required to find a solution. However, in
the general nonlinear case, several iterations of Newton’s step are required.
A common approach to guarantee the convergence of the overall iterative
procedure is to introduce a line search and a merit function [22]. Given a
guess at iteration $i$ and a direction $p$, the next guess is defined by

$\displaystyle\begin{pmatrix}x^{i+1}_{0:T}\\\
u^{i+1}_{t:T-1}\end{pmatrix}=\begin{pmatrix}x^{i}_{0:T}\\\
u^{i}_{t:T-1}\end{pmatrix}+\alpha_{i}\begin{pmatrix}p_{x_{0:T}}\\\
p_{u_{t:T-1}}\end{pmatrix},$ (15)

where the step length $\alpha_{i}$ is chosen in order to decrease the merit
function. As advocated by Nocedal et al. [22], a Newton step provides a
descent direction for the merit function

$f_{\mathcal{M}}(x_{0:T},u_{t:T-1})=\frac{1}{2}\sum_{j=0}^{T}\left\|\frac{\partial
J}{\partial
x_{j}}\right\|^{2}+\frac{1}{2}\sum_{j=t}^{T-1}\left\|\frac{\partial
J}{\partial u_{j}}\right\|^{2}.$ (16)

This merit function can be derived analytically. The exact derivations of each
gradient are provided in the supplementary material [23]. By construction, the
expected decrease of the direction $p$ derived by Algorithm 1 is $p^{T}\nabla
f_{\mathcal{M}}=-||\nabla J||_{2}^{2}$ and the following convergence
guarantees hold.

######  Proposition 4.

Assuming that the norm of the inverse of the Hessian of $J$ is bounded and
that the step length $\alpha_{i}$ satisfies the Wolfe conditions for the merit
function (16), the sequence $(x^{i}_{0:T},u^{i}_{t:T-1})_{i}$ defined by the
update rule (15) with steps from Algorithm 1 is globally convergent to a
stationary point of (3). Furthermore, when the iterate is sufficiently close
to the solution, the sequence has quadratic convergence.

###### Proof.

This proposition is a direct consequence from the fact that Algorithm 1 yields
a Newton step. A detailed proof of the convergence guarantees of Newton method
is provided in [22]. ∎

One may ask under which conditions the Hessian of $J$ is non-degenerate.
Intuitively, large values of $\mu$ can make the problem ill-defined. Indeed,
if the opposing player can choose large disturbances, then the controller
might not be able to compensate. In the linear quadratic case, Whittle [8]
studied this maximum value for $\mu$ that makes the problem ill-defined.
Although, we do not study this limit value in the nonlinear case, we note that
analogously to the linear quadratic case, the algorithm is well defined if
$\Gamma_{k}$ and ${P_{k}^{-1}-\mu\bar{\ell}_{k}^{xx}}$ are positive-definite
which is the case when $\mu$ is small enough.

###  III-E About the cooperative case

So far, only the case, $\mu>0$ has been considered, but, the case $\mu<0$ is
also well defined. Indeed, the search for a stationary point of $J$ can be
done for an arbitrary sign of $\mu$. However, such a stationary point would
now be a way to find a local minimum of $J$ with respect to the variables
${(x_{0:T},u_{t:T-1})}$. Interestingly, this scenario can be interpreted as a
cooperative scenario between the controller and the opposing player. In fact,
the disturbances can be seen a second controller minimizing the cost, $\ell$,
and maximizing the likelihood of the disturbances. Clearly, this change of
sign does not affect the derivation of the stagewise Newton method. However,
it can be noted that in that case, one can directly use the cost $J$ as a
merit function.

##  IV Experiments

A Python implementation of the proposed method is available online [23]. It is
based on the Crocoddyl software [6], a state of the art (risk-neutral) DDP
solver which provides analytical derivatives of robot dynamics. In this
section, two numerical examples illustrate the proposed method.

###  IV-A Planar quadrotor

We study a quadrotor moving in a plane aiming to reach the position $(p_{x}\ \
p_{y})=(2\ \ 0)$ starting at the origin without initial velocity. The state is
${x=(p_{x}\ \ p_{y}\ \ \theta\ \ \dot{p}_{x}\ \ \dot{p}_{y}\ \
\dot{\theta})^{T}}$ where $\theta$ is the orientation of the quadrotor. The
system dynamics is

$\displaystyle m\ddot{p}_{x}$ $\displaystyle=-(u_{1}+u_{2})\sin(\theta),$
$\displaystyle m\ddot{p}_{y}$ $\displaystyle=(u_{1}+u_{2})\cos(\theta)-mg,$
$\displaystyle J\ddot{\theta}$ $\displaystyle=r(u_{1}-u_{2}),$ (17)

where the control input $u=(u_{1}\ \ u_{2})^{T}\in\mathbb{R}^{2}$ represents
the force at each rotor and $g$ is gravitational acceleration. An exponential
cost models the presence of an obstacle

$\displaystyle\ell_{k}(x_{k},u_{k})$
$\displaystyle=0.3\exp\left(-10({p_{x}}_{k}-1)^{2}-0.5({p_{y}}_{k}+0.1)^{2}\right)$
$\displaystyle+0.005\|u_{k}-\bar{u}\|_{2}^{2}+0.05(x_{k}-x_{\star})^{T}L(x_{k}-x_{\star})$
$\displaystyle\ell_{T}(x_{T})$
$\displaystyle=(x_{T}-x_{\star})^{T}L(x_{T}-x_{\star}),$ (18)

where $\bar{u}=\frac{mg}{2}(1\ \ 1)^{T}$, $x_{\star}=(2\ \ 0\ \ 0\ \ 0\ \ 0\ \
0)^{T}$, and $L=\operatorname{diag}(100\ \ 100\ \ 100\ \ 1\ \ 1\ \ 1)$. Only
the position, $p_{x}$, $p_{y}$, and orientation $\theta$ are part of the
measurements. The integration and discretization of the model is done with
Runge Kutta 4 with a time step of $0.05$ and the total horizon, $T$, is $60$.
Furthermore, ${P_{0}=Q=10^{-5}I_{6}}$, ${R=10^{-4}\operatorname{diag}(1\ \ 1\
\ 0.01)}$ and $\hat{x}_{0}$ is the origin. A backtracking line-search is used
– a step is accepted if ${f_{\mathcal{M}}^{i+1}\leq
f_{\mathcal{M}}^{i}+\alpha_{i}cp_{i}^{T}f_{\mathcal{M}}^{i}}$ where
${f_{\mathcal{M}}^{i}=f_{\mathcal{M}}(x^{i}_{0:T},u^{i}_{t:T-1})}$ with
${c=\frac{1}{4}}$. The iterative process is stopped when the decrease in the
merit function is lower than $10^{-12}$.

Figure 1 illustrates the solution obtained by the solver for different values
of $\mu$. The neutral controller, the limit when $\mu$ tends to zero, aiming
at minimizing the cost without accounting for disturbances, is solved using
DDP. Interestingly, when $\mu$ is positive (the non-cooperative case), the
opposing player chooses disturbances that will push the quadrotor towards to
the obstacle and when $\mu$ is negative (the cooperative case), the opposing
player chooses disturbances that will push the quadrotor away form the
obstacle. For this experiment, we found that the value of $\mu$ for which
$\Gamma_{k}$ and ${P_{k}^{-1}-\mu\bar{\ell}^{xx}_{k}}$ are no longer positive-
definite matrices was around $20$.

Figure 1: Initial plan for different values of $\mu$. The larger $\mu$ the
more the controller plans to be pushed against the obstacle. Figure 2: Average
trajectory. Compared the neutral controller, the dynamic game controller
($\mu=6$) exhibits a risk sensitive behavior as it remains further from the
high cost area representing the obstacle.  Figure 3: Average control
trajectories. Compared the neutral controller, the dynamic game controller
($\mu=6$) has a larger standard deviation.

Next, we illustrate the risk-sensitive behavior of the resulting controller in
closed loop in a simulation with disturbances following Gaussian
distributions: ${x_{0}\sim\mathcal{N}(\hat{x}_{0},P_{0})}$,
${w_{k}\sim\mathcal{N}(0,Q)}$ and ${\gamma_{k}\sim\mathcal{N}(0,R)}$. We set
$\mu=6$ and at each time step of the simulation, Algorithm 1 is run and
$u_{t}$ is applied to the system. Additionally, we compare to the neutral
controller – an iterative Kalman smoother is used for the filtering and the
controller uses DDP with the last state estimate from the smoother as an
initial condition. Figure 2 depicts the average and standard deviation over
one thousand simulation. We can see that, in this MPC scheme, the dynamic game
controller maintains a larger distance from the obstacle, resulting in safer
behavior. Figure 3 shows the distribution of control trajectories.
Interestingly, the dynamic game controller has a larger standard deviation.

###  IV-B An industrial robot

In this example, we consider the 7-DoF torque-controlled KUKA LWR iiwa R820
14. The dynamics of the robot are provided by Pinocchio [24]. The
14-dimensional state is composed of the joint positions and velocities. We
consider the following prior on the initial condition
$\hat{x}_{0}=\begin{pmatrix}0.1&0.7&0.&0.7&-0.5&1.5&0.\end{pmatrix}^{T}$. The
control input is a 7-dimensional vector of the torque applied on each joint.
The goal is to move the end effector to a desired position,
$p_{\text{target}}=\begin{pmatrix}-0.4&0.3&0.7\end{pmatrix}$ with the
following cost:

$\displaystyle\ell_{k}(x_{k},u_{k})$
$\displaystyle=10^{-3}\|x_{k}-\hat{x}_{0}\|_{2}^{2}+10^{-6}\|u_{k}-\bar{u}(x_{k})\|_{2}^{2}$
$\displaystyle+10^{-1}\|p_{\text{target}}-\bar{p}(x_{k})\|_{2}^{2}$
$\displaystyle\ell_{T}(x_{T})$
$\displaystyle=\|p_{\text{target}}-\bar{p}(x_{k})\|_{2}^{2}+10^{-3}\|x_{k}-\hat{x}_{0}\|_{2}^{2},$
(19)

where $\bar{u}(x_{k})$ is the gravity compensation and $\bar{p}(x_{k})$ the
position of the end-effector. For the measurement model, we assume that only
the joints position are measured. The initial control inputs $u_{0},\dots
u_{t-1}$ are generated with DDP. Given those $t$ initial control inputs, $t$
measurements are generated according to an undisturbed trajectory. More
precisely, for $1\leq k\leq t$, the observations are defined by
${y_{k}=h_{k}(f^{(k)}(\hat{x}_{0},u_{1:k-1}))}$ where
$f^{(k)}(\hat{x}_{0},u_{1:k-1})$ denotes the state value at the $k^{th}$ step
integrated from the initial guess $\hat{x}_{0}$. The horizon, $T$, is $100$
and $t$ is $5$ while the sensitivity parameter is set to $\mu=1.2$. Here
$P_{0}=Q=0.01\times I_{14}$ and $R=0.5\times I_{7}$. For this experiment, the
second order derivatives of the dynamics were approximated to be null as the
former are not provided by most standard rigid body dynamics libraries such as
Pinocchio [24]. Note that this is a common practice for state of the art
optimal control algorithms in robotics [6]. Our solver converges nevertheless,
suggesting that this approximation might be used to scale to a large number
degrees of freedom for real-time computations (e.g. MPC). In Figure 4, we plot
the solution of the dynamics game solver compared to DDP in the end effector
space, the shaded grey area represents the estimation part of the solution. We
can see that the dynamic game controller plans that the disturbances will slow
down the reaching task.

Figure 4: End-effector position $p^{EE}$ vs time. The dashed lines represent
the target.

##  V Conclusion

We introduced an iterative solver to find local Nash equilibrium of dynamic
game with imperfect state measurements. The proposed algorithm is proven to be
equivalent to Newton’s method and benefits from its convergence properties
while scaling linearly with the time horizon.

## References

  * [1] D. Q. Mayne, “Model predictive control: Recent developments and future promise,” _Automatica_ , vol. 50, pp. 2967–2986, 2014.
  * [2] M. R. James, J. S. Baras, and R. J. Elliott, “Risk-sensitive control and dynamic games for partially observed discrete-time nonlinear systems,” _IEEE transactions on automatic control_ , vol. 39, pp. 780–792, 1994.
  * [3] T. Başar and P. Bernhard, “$\mbox{H}^{\infty}$-0ptimal control and related minimax design problems: A dynamic game approach,” _IEEE Trans. Autom. Control._ , vol. 41, 1996.
  * [4] D. A. Copp and J. P. Hespanha, “Nonlinear output-feedback model predictive control with moving horizon estimation,” in _53rd IEEE Conference on Decision and Control_. IEEE, 2014, pp. 3511–3517.
  * [5] ——, “Simultaneous nonlinear model predictive control and state estimation,” _Automatica_ , vol. 77, pp. 143–154, 2017.
  * [6] C. Mastalli, R. Budhiraja, W. Merkt, G. Saurel, B. Hammoud, M. Naveau, J. Carpentier, L. Righetti, S. Vijayakumar, and N. Mansard, “Crocoddyl: An Efficient and Versatile Framework for Multi-Contact Optimal Control,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2020.
  * [7] D. Jacobson, “Optimal stochastic linear systems with exponential performance criteria and their relation to deterministic differential games,” _IEEE Transactions on Automatic control_ , vol. 18, pp. 124–131, 1973.
  * [8] P. Whittle, “Risk-sensitive linear/quadratic/gaussian control,” _Advances in Applied Probability_ , vol. 13, pp. 764–777, 1981.
  * [9] B. Hammoud, A. Jordana, and L. Righetti, “irisc: Iterative risk sensitive control for nonlinear systems with imperfect observations,” _arXiv preprint arXiv:2110.06700_ , 2021.
  * [10] M. C. Campi and M. R. James, “Nonlinear discrete-time risk-sensitive optimal control,” _International Journal of Robust and Nonlinear Control_ , vol. 6, pp. 1–19, 1996.
  * [11] T. Başar, “Robust designs through risk sensitivity: An overview,” _Journal of Systems Science and Complexity_ , vol. 34, pp. 1634–1665, 2021\.
  * [12] P. Bernhard, “Minimax versus stochastic partial information control,” in _Proceedings of 1994 33rd IEEE Conference on Decision and Control_ , vol. 3, 1994, pp. 2572–2577 vol.3.
  * [13] H. Cox, “On the estimation of state variables and parameters for noisy dynamic systems,” _IEEE Transactions on automatic control_ , vol. 9, pp. 5–12, 1964\.
  * [14] H. E. Rauch, F. Tung, and C. T. Striebel, “Maximum likelihood estimates of linear dynamic systems,” _AIAA journal_ , vol. 3, pp. 1445–1450, 1965.
  * [15] B. M. Bell, “The iterated kalman smoother as a gauss–newton method,” _SIAM Journal on Optimization_ , vol. 4, pp. 626–636, 1994.
  * [16] D. Murray and S. Yakowitz, “Differential dynamic programming and newton’s method for discrete optimal control problems,” _Journal of Optimization Theory and Applications_ , vol. 43, pp. 395–414, 1984.
  * [17] J. C. Dunn and D. P. Bertsekas, “Efficient dynamic programming implementations of newton’s method for unconstrained optimal control problems,” _Journal of Optimization Theory and Applications_ , vol. 63, pp. 23–38, 1989\.
  * [18] H. Mukai, A. Tanikawa, I. Tunay, I. Ozcan, I. Katz, H. Schattler, P. Rinaldi, G. Wang, L. Yang, and Y. Sawada, “Game-theoretic linear quadratic method for air mission control,” in _39th IEEE Conference on Decision and Control_ , vol. 3. IEEE, 2000, pp. 2574–2580.
  * [19] J. Morimoto, G. Zeglin, and C. Atkeson, “Minimax differential dynamic programming: application to a biped walking robot,” in _Proceedings 2003 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2003)_ , vol. 2, 2003, pp. 1927–1932.
  * [20] B. Di and A. Lamperski, “Newton’s method, bellman recursion and differential dynamic programming for unconstrained nonlinear dynamic games,” _Dynamic Games and Applications_ , pp. 1–49, 2021.
  * [21] T. Başar and G. J. Olsder, _Dynamic noncooperative game theory_. SIAM, 1998.
  * [22] J. Nocedal and S. J. Wright, _Numerical optimization_. Springer, 1999.
  * [23] A. Jordana, B. Hammoud, J. Carpentier, and L. Righetti,
https://github.com/machines-in-motion/dynamic˙game˙optimizer.

  * [24] J. Carpentier, F. Valenza, N. Mansard, _et al._ , “Pinocchio: fast forward and inverse dynamics for poly-articulated systems,” https://stack-of-tasks.github.io/pinocchio, 2015–2021.