* Mandelbaum et al. (2013) Mandelbaum R., Slosar A., Baldauf T., Seljak U., Hirata C. M., Nakajima R., Reyes R., Smith R. E., 2013, MNRAS, 432, 1544
  * Marques et al. (2020) Marques G. A., Liu J., Huffenberger K. M., Colin Hill J., 2020, ApJ, 904, 182
  * McEwen et al. (2016) McEwen J. E., Fang X., Hirata C. M., Blazek J. A., 2016, J. Cosmology Astropart. Phys, 2016, 015
  * Miyatake et al. (2022) Miyatake H., et al., 2022, Phys. Rev. D, 106, 083520
  * Miyatake et al. (2023) Miyatake H., et al., 2023, Phys. Rev. D, 108, 123517
  * Myles et al. (2021) Myles J., Alarcon A., et al., 2021, MNRAS, 505, 4249
  * Namikawa et al. (2019) Namikawa T., et al., 2019, ApJ, 882, 62
  * Omori et al. (2023) Omori Y., et al., 2023, Phys. Rev. D, 107, 023529
  * Pandey et al. (2020) Pandey S., et al., 2020, Phys. Rev. D, 102, 123522
  * Pandey et al. (2022) Pandey S., et al., 2022, Phys. Rev. D, 106, 043520
  * Porredon et al. (2021) Porredon A., et al., 2021, Phys. Rev. D, 103, 043503
  * Porredon et al. (2022) Porredon A., et al., 2022, Phys. Rev. D, 106, 103530
  * Prat et al. (2018) Prat J., Sánchez C., et al., 2018, Phys. Rev. D, 98, 042005
  * Prat et al. (2022) Prat J., et al., 2022, Phys. Rev. D, 105, 083528
  * Prat et al. (2023a) Prat J., Zuntz J., et al., 2023a, The Open Journal of Astrophysics, 6, 13
  * Prat et al. (2023b) Prat J., Zacharegkas G., et al., 2023b, MNRAS, 522, 412
  * Raveri & Hu (2019) Raveri M., Hu W., 2019, Phys. Rev. D, 99, 043506
  * Robertson et al. (2021) Robertson N. C., et al., 2021, A&A, 649, A146
  * Rodríguez-Monroy et al. (2022) Rodríguez-Monroy M., et al., 2022, MNRAS, 511, 2665
  * Samuroff et al. (2017) Samuroff S., Troxel M. A., Bridle S. L., Zuntz J., MacCrann N., Krause E., Eifler T., Kirk D., 2017, MNRAS, 465, L20
  * Samuroff et al. (2019) Samuroff S., et al., 2019, MNRAS, 489, 5453
  * Samuroff et al. (2023) Samuroff S., et al., 2023, MNRAS, 524, 2195
  * Sánchez et al. (2022) Sánchez C., Prat J., et al., 2022, Phys. Rev. D, 105, 083529
  * Sánchez et al. (2023) Sánchez C., Alarcon A., et al., 2023, MNRAS, 525, 3896
  * Schneider et al. (2000) Schneider P., King L., Erben T., 2000, A&A, 353, 41
  * Secco et al. (2022) Secco L. F., Samuroff S., et al., 2022, Phys. Rev. D, 105, 023515
  * Serrano et al. (2023) Serrano S., et al., 2023, MNRAS, 523, 3287
  * Sheldon et al. (2004) Sheldon E. S., et al., 2004, AJ, 127, 2544
  * Singh et al. (2015) Singh S., Mandelbaum R., More S., 2015, MNRAS, 450, 2195
  * Singh et al. (2020) Singh S., Mandelbaum R., Seljak U., Rodríguez-Torres S., Slosar A., 2020, MNRAS, 491, 51
  * Sugiyama et al. (2023) Sugiyama S., et al., 2023, Phys. Rev. D, 108, 123521
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152
  * Tröster et al. (2022) Tröster T., et al., 2022, A&A, 660, A27
  * Troxel et al. (2018) Troxel M. A., et al., 2018, Phys. Rev. D, 98, 043528
  * van Uitert et al. (2018) van Uitert E., et al., 2018, MNRAS, 476, 4662
  * Yoo & Seljak (2012) Yoo J., Seljak U., 2012, Phys. Rev. D, 86, 083504
  * Yoon et al. (2019) Yoon M., Jee M. J., Tyson J. A., Schmidt S., Wittman D., Choi A., 2019, ApJ, 870, 111
  * Zacharegkas et al. (2022) Zacharegkas G., et al., 2022, MNRAS, 509, 3119
  * Zhou et al. (2021) Zhou R., et al., 2021, MNRAS, 501, 3309
  * Zuntz et al. (2015) Zuntz J., et al., 2015, Astronomy and Computing, 12, 45

##  Appendix A Appendix A. Posteriors on redshift, galaxy bias and cosmology
parameters

Figure 13: Posteriors on cosmology, redshift and galaxy bias parameters for
the example IA scenario discussed in Section 4.1. As in that section, we
choose to focus on redMaGiC lenses here, but a similar picture is seen with
MagLim. The purple and black contours show cosmic shear and galaxy-galaxy
lensing + clustering analyses respectively. Although there are slight
($<1\sigma$) shifts in the lowest source redshift bin (labelled $\Delta
z^{1}_{\rm s}$ here), in general the source and lens redshift parameters and
galaxy bias are stable and close to their input values (indicated by the
dotted lines). The intrinsic alignment amplitude $A_{1}$, on the other hand,
is seen to shift significantly in the $2\times 2$pt case. Figure 14:
Posteriors on the full cosmology space from our mock cosmic shear only and
galaxy-galaxy lensing plus clustering analyses. Overlain in blue is the prior
(note this is a set of $\sim 10,000$ samples drawn from a $>15$ dimensional
prior, and so there is some noise). The fine dotted lines in each panel show
the input values used to generate the data vector. Although the cosmic shear
only analysis (purple) shows some bias in this space, with $h$,
$n_{\mathrm{s}}$ and $\Omega_{\mathrm{b}}$ pushing up against the prior edges,
this is much less true for either $2\times 2$pt posteriors (black shaded and
dotted), which are largely unbiased.

This appendix presents the posteriors in the parameter space of cosmology, IAs
and redshift error parameters. Figure 13 shows cosmic shear only (purple) and
galaxy-galaxy lensing + galaxy clustering (black) analyses of our example data
vector from Section 4.1. The input parameters are indicated by the dashed
lines. Although the IA parameters do shift, particularly in the $2\times 2$pt
case, the redshift parameters are noticeably stable. In neither analysis are
either the lens or the source redshift nuisance parameters shifting to absorb
the IA error present in the data. This is most likely at least in part due to
the priors – these are Gaussian and relatively informative, which restricts
the ability of these parameters to vary too widely (see Table 1). Note that we
are using redMaGiC here, but a similar plot using MagLim lenses paints the
same overall picture. Even with the extra width parameters associated with
MagLim, we do not see a significant amount of IA bias being absorbed into the
redshift distributions.

Figure 14 shows the wider cosmology space. Overlain in blue is the prior. In
the main body of the paper we focus on $S_{8}$ and $\Omega_{\mathrm{m}}$,
since these are the most constrained parameters. It can also, however, be
interesting to look at the behaviour of the other, less constrained,
parameters. Again, we do not see an obvious explanation for the relative lack
of $S_{8}$ bias in the $2\times 2$pt case. The neutrino mass parameter
$\Omega_{\nu}h^{2}$ is almost completely unaffected, being strongly prior
dominated. In general, $2\times 2$pt (black filled and dashed) is in fact
slightly more stable throughout the parameter space than cosmic shear
(purple), which is shifted towards the prior edge in
$\Omega_{\mathrm{b}},n_{\mathrm{s}}$ and $h$.

The overall picture here is fairly simple – $2\times 2$pt appears to be less
biased in general than $1\times 2$pt, throughout the whole cosmological
parameter space, and not just the projected $S_{8}-\Omega_{\mathrm{m}}$ plane.
Interestingly the IA modelling error appears to be absorbed primarily by the
IA parameters, without major leakage into other parts of parameter space.

##  Appendix B Appendix B. Comparing lens samples

In this appendix we consider in more detail the differences between our
redMaGiC and MagLim lens samples, and the impact those differences have on our
results.

Since the choice of lens sample has some effect on the level of cosmological
bias in a NLA-based $2\times 2$pt analysis (see Figure 4), it is worth trying
to unwrap the differences. That is, going from redMaGiC to MagLim groups
together a few different changes, and it is useful to consider them
separately. The major differences are (a) the number of lens bins (5 for
redMaGiC and either 4 or 6 for MagLim); (b) the shape of the estimated $n(z)$
themselves (c) the total effective number density $n_{\rm eff}$ of our MagLim
lens sample is greater than in redMaGiC ($n_{\rm eff}=0.72$ galaxies/square
arcmin and 0.17); (d) the uncertainties on the redshift estimates are larger
for MagLim (due to the inclusion of fainter galaxies), and so the width of the
priors on redshift nuisance parameters is slightly larger. We consider each of
these in turn below.

Beginning with (a), the number of lens bins clearly has an impact on
constraining power. This can be seen in Figure 4 – the pink contour is
noticeably smaller than the purple. Including the upper two MagLim bins
clearly helps to constrain cosmology. It does not, however, shift the contour
back to sit on top of the black (redMaGiC) one. Although the test is not
perfect (since different bins have different redshift priors), it does imply
that the choice of binning alone is not the main factor in the redMaGiC-MagLim
differences.

Likewise we do not expect (b) to be the primary cause of differences in IA
sensitivity. Although in principle differences in the shape and mean redshift
of the lens distributions does alter the lens-source overlap, and thus the
sensitivity to IAs, in practice the DES Y3 MagLim and redMaGiC $n(z)$s are not
hugely different. This can be seen qualitatively by comparing the lower two
panels of Figure 1. A more precise comparison is given in Figure 15, where we
show the integrated source-lens overlap, $X=\int n_{\rm source}(z)n_{\rm
lens}(z)\mathrm{d}z$, which is the effective redshift kernel entering the
galaxy-intrinsic ($g$I) contribution to $\gamma_{t}$. Again the differences
are relatively small (compare the upper and lower panel; note that the y-axes
are the same). On the other hand, the combination of (c) and (d) _is_ thought
to play a more significant part in the observed differences. These factors
taken together are already known to have a significant impact on the $2\times
2$pt cosmological constraining power (see Porredon et al. 2021 Figure 10).
That is, the increased density outweighs the slightly broader redshift priors,
resulting in a tighter constraint in the $S_{8}-\Omega_{\mathrm{m}}$ plane.
This added density reduces the statistical uncertainty, and thus reduces the
tolerance of small biases. Considering Figure 4, however, the difference
between the black and purple contour is not simply a contraction about a
shifted mean; a reduction in the uncertainty on $S_{8}$ by $10\%$ would not
ultimately make much difference given the relatively accurate centring. The
differences seen in Figure 4 are thought to primarily arise from the broader
$\Delta z$ and $\sigma z$ priors. In several bins we go from the $n(z)$ width
being fixed entirely to being varied with a Gaussian prior (see the lower part
of Table 1). These width parameters are known to be degenerate with galaxy
bias, and so with the overall amplitude of the $2\times 2$pt data. By
modifying the lens-source overlap, these parameters can effectively control
how a given IA signal enters the data. This creates regions of parameter space
where $S_{8}$ and IA error are degenerate, and so reduces the ability of
different parts of the $\gamma_{t}$ data to self-calibrate such modelling
error.

Figure 15: A quantification of the lens-source overlap for each possible bin
pair. The quantity on the vertical axis is the total area enclosed by the
product of the lens and source redshift distributions; this overlap term
effectively modulates the size of the IA signal entering our $\gamma_{t}$
predictions (see the equations in Section 2). Each coloured band shows a
particular lens bin. The upper and lower panels show the same calculation for
redMaGiC and MagLim lens samples respectively. The overall pattern is seen to
be roughly the same between the two lens samples. Pairs involving the upper
two MagLim lens bins are represented by open circles, as these are not
included in the fiducial DES Y3 analysis. Figure 16: The fractional
contributions of the $\delta_{g}$I and $\mu$I (galaxy-intrinsic and
magnification-intrinsic) IA terms to our galaxy-galaxy lensing data. Both are
shown relative to the cosmological signal (the sum of galaxy-shear and
magnification-shear contributions $\delta_{g}$G and $\mu$G). As in Figure 9,
shaded grey bands represent scale cuts. The numbers shown in each panel
specify a particular bin pair $(l,s)$. The solid lines show the input IA
signal for our mock $\gamma_{t}$ data, and dashed lines show the best fit
prediction from an NLA $2\times 2$pt analysis of that data. Note we are
showing the MagLim lens configuration here. See Figure 10 for the equivalent
of this plot using the redMaGiC lens sample.

##  Appendix C Appendix C. Building an LSST Y1 like simulated analysis

In this appendix we describe in more detail the details of our LSST Y1 like
mock analyses. These are used in Section 4 to test the robustness of our
results in an at least semi-realistic Stage IV like setup. The idea is not to
predict every aspect of the future Rubin analysis correctly, but more to test
if our findings are still applicable after a significant gain in constraining
power.

Bin | $n_{\rm eff}$ | Bias $b_{1}$ | $C$
---|---|---|---
1 | 3.60 | 1.42 | 0.43
2 | 3.60 | 1.66 | 0.30
3 | 3.60 | 1.70 | 1.75
4 | 3.60 | 1.62 | 1.94
5 | 3.60 | 1.78 | 1.56
Table 2: Lens sample parameters for our LSST Y1 like mock data. Shown (left to
right) for each redshift bin are the effective number density, linear galaxy
bias and magnification coefficient.

For our LSST Y1 source sample, we follow Prat, Zuntz et al. (2023a) and define
5 equal number redshift bins. We assume the true underlying distribution has
the form $p(z)\propto z^{2}\times\mathrm{exp}[-(z/z_{0})^{\alpha}]$, where
$\alpha=0.78$ and $z_{0}=0.13$ (these numbers are taken from the DESC SRD -
i.e. DESC Collaboration 2018’s Sec D2.1 and Appendix F2). A lower limit of
$z>0.2$ is imposed, sharp bin edges are defined, and the resulting
distributions are convolved with a Gaussian redshift error $\sigma_{z}=0.05$.
This results in the source $n(z)$s shown in the upper panel of Figure 17
(which approximately match Figure 3 of Prat, Zuntz et al. 2023a). We follow
DESC Collaboration 2018 and Fang et al. (2020) in assuming a total number
density of 11.2 galaxies per square arcminute, which divided across 5 bins,
gives $2.24$ per bin. We also assume an ellipticity dispersion of
$\sigma_{e}=0.26$ in all bins.

Figure 17: Redshift distributions for our mock LSST Y1 like samples. Note that
these follow the prescription of DESC Collaboration (2018). See Appendix C for
discussion.

For our lens sample, we assume five bins and again use an analytic form for
the redshift true distribution ($\alpha=0.94$, $z_{0}=0.26$; see DESC
Collaboration 2018). Following the predictions in the SRD, we truncate this
distribution at $0.2<z<1.2$ and apply a Gaussian error of
$\sigma_{z}=0.03\times(1+z)$. The resulting lens $n(z)$ in five bins are shown
in the lower panel of Figure 17. We keep the Y3 MagLim bias and magnification
values as our inputs, as shown in Table 2. For the lens sample number density,
we adopt the DESC Collaboration (2018) prediction for a Y1 photometric sample
of 18 galaxies per square arcminute. Divided equally between our 5 bins, this
gives us 3.60 per square arcmin.

Note that in both cases, but particularly the source sample, these estimates
will likely be quite different from the distributions eventually estimated
from real LSST Y1 data. For one thing, the idealised selection function of
DESC Collaboration (2018) will likely become more complicated, which could
affect the depth as well as the shape of the distributions. Secondly, modern
photometric redshift algorithms, designed to robustly characterise the
ensemble $n(z)$s for a particular population (see e.g. Myles, Alarcon et al.
2021, Campos et al. 2024) typically result in non-Gaussian distributions with
asymmetric tails to high/low redshift. Underestimating the tails is not
necessarily an issue for predicting the shear signal alone, but it does matter
for the IA contribution (sensitive to the shape of the lens and source
$n(z)$s). We should, then, consider our forecasts as a rough exercise to see
how things change with greater depth and density rather than a robust
prediction for LSST Y1.

Finally, we assume a joint lens-source mask with an area of 12,300 square
degrees. The joint covariance matrix of our LSST Y1 like shear, galaxy-galaxy
lensing and clustering data is calculated using CosmoCov. This computes both
Gaussian shape noise and cosmic variance terms, as well as super sample and
connected non-Gaussian contributions. Note this is the same methodology as was
used for DES Collaboration (2022), albeit with updated parameters.

For the analysis itself, we use the DES Y3 pipeline described in Section 2,
with some differences. While source sample and cosmology priors are kept the
same as in Table 1, we update the lens sample priors to create a mock sample
with redMaGiC-comparable photometry. That is, we keep the width parameters
fixed and adopt Gaussian priors on the $\Delta z$ parameters of widths 0.004,
0.003, 0.003, 0.005 and 0.01 respectively in our five lens bins.

For scale cuts, we apply the DES Y3 method of Krause et al. (2021). That is,
we generate mock LSST Y1 like data with and without baryonic contamination
from the OWLS-AGN simulations (van Daalen et al., 2011). Scales are
iteratively removed from the joint $\xi_{\pm}$ data vector until the total
$\chi^{2}$ between the two data vectors $(\mathbf{D}^{\rm
OWLS}-\mathbf{D}^{\rm fid})\mathbf{C}^{-1}(\mathbf{D}^{\rm
OWLS}-\mathbf{D}^{\rm fid})$ is below some threshold value. For the purposes
of this work, we choose $\chi^{2}=2.5$ – this is slightly less strict than the
value eventually adopted by the DES collaboration (DES Collaboration, 2022),
where the fiducial cuts were defined by $\chi^{2}=0.4$. There are two
competing factors to account for here. The first, and perhaps more easy to
predict, is that LSST Y1 will simply have more constraining power than
previous surveys. This tightening of statistical uncertainties will allow less
tolerance for a given degree of systematic error, and so will ultimately
require more stringent scale cuts for any given model. The second is that
there will inevitably be advances in the field of baryon modelling in the
coming years. Even if it is not trusted on _all_ scales, it seems likely that
Rubin will incorporate some sort of baryonic $P(k)$ model, and so access
slightly smaller scales than they would if they were relying on a dark matter
only model. It was this logic that led to our choice of 2.5. We should bear in
mind that this is not designed as a robust prediction, but at least roughly
accounts for these competing factors. For the galaxy-galaxy lensing and
clustering parts of the data vector, we again follow Krause et al. (2021).
That is, we convert limits in comoving space at $6$ Mpc$/h$ (for $\gamma_{t}$)
and $8$ Mpc$/h$ (for $w$) into angular scales using the mean redshift of each
lens bin. We run two types of cosmological analyses in this LSST Y1 setup. The
first is $\Lambda$CDM, which uses the same priors and parameter choices as
shown in Table 1. The second is a $w$CDM analysis, which is the same, but with
an extra equation of state parameter $w$, which is varied with a flat prior
over the range $\mathrm{U}[-3,-0.33]$.

One consequence of the extra constraining power is that PolyChord takes
considerably longer to converge. For this reason, for the LSST like chains we
swap to Nautilus. The basic idea of this new sampler is to use neural networks
to help choose boundaries in parameter space most efficiently, and so minimise
the total number of likelihood calls needed to estimate the posterior. The
algorithm comes packaged with CosmoSIS v3, and can be called in the same way
as other samplers. We validate the convergence and accuracy of our LSST Y1
Nautilus chains by re-running one of them using a more conventional MCMC
approach – emcee (Foreman-Mackey et al., 2013). Specifically we re-run our
less extreme IA scenario $1\times 2$pt analysis (the purple contour in the
upper right panel of Figure 11). After allowing our emcee chain to run for
just over 1M samples, we cut the first quarter for burn-in. Plotting the
surviving 750,000, we find good agreement with the shorter Nautilus chain (the
marginalised mean and standard deviation of $S_{8}$ are consistent to $<1\%$).
Note that we do not use the Bayesian evidence values from Nautilus in this
work, so we do not compare or try to validate those here.

##  Appendix D Appendix D. The impact of shear ratios

Figure 18: The impact of shear ratios. We show our fiducial $1\times 2$pt and
$2\times 2$pt posteriors as shaded purple/black contours. The dashed contours
are the same analyses, but including an additional likelihood from shear
ratios on scales $2-6$ Mpc$/h$. Note that the simulated shear ratio data
vector contains the same input IA scenario as the large scale $\xi_{\pm}$ and
$\gamma_{t}$ data. As in previous plots, the input parameters are marked with
a star.

In this appendix we briefly consider the impact of shear ratios (SR) on our
results. SR are the ratios of galaxy-galaxy lensing measurements using the
same lens bin. In theory the power spectra should (at least partially) cancel
away, leaving ratios of angular diameter distances (see Jain & Taylor 2003,
Sánchez, Prat et al. 2022). For this reason they were proposed as a way to
calibrate redshift error (Prat et al., 2018). In practice, they can also
contain a significant amount of information about intrinsic alignments (e.g.
Blazek et al. 2012). The fiducial DES Y3 cosmic shear, $2\times 2$pt and
$3\times 2$pt analyses (Secco, Samuroff et al. 2022; Amon et al. 2022;
Porredon et al. 2022; DES Collaboration 2022) all included SR. We chose not to
for the main sections of this paper because doing so mixes information coming
from $\gamma_{t}$ and $\xi_{\pm}$, and so complicates our interpretation of
differences between the probes. Since SR potentially carry information about
IAs, however, it is interesting to briefly consider how our main results
change as they are added.

To avoid significant correlation with the large scale $\gamma_{t}$ and
$\xi_{\pm}$ data, the DES Y3 setup includes SR measurements on scales of $2-6$
Mpc$/h$, where the uncertainties are shape noise dominated. Note that this is
outside the one-halo regime, but roughly where higher-order TATT terms are
expected to become significant (see Blazek et al. 2019). The covariance of the
SR is propagated from the analytic covariance matrix used for the large scale
inference. The SR data are averaged across angular scales, giving one data
point per lens-source-source triplet. In total we have 9 data points from SR,
constituting the first three lens bins, each paired with source bins (1,4),
(2,4) and (3,4) (Sánchez, Prat et al., 2022). The Y3 pipeline incorporates
them as an additional likelihood, such that
$\mathrm{ln}\mathcal{L}=\mathrm{ln}\mathcal{L}_{\rm
2pt}+\mathrm{ln}\mathcal{L}_{\rm SR}$ (see e.g. Secco, Samuroff et al. 2022
Section IV-A).

To test how this extra information would interact with IA error, we set up a
simulated test along the lines of those presented in Section 4. Now, however,
we also include a simulated SR data vector containing our fiducial TATT
scenario (see Section 4.1). We rerun our chains in this setup, modelling both
the large scale two-point functions and the SR using NLA. We should note that
this is a slightly contrived setup to help us understand how SR affect the
dynamics of self calibration. The fiducial DES Y3 chains used TATT, which is
valid on scales well below NLA.

Our results are shown in Figure 18. As we can see, the cosmic shear posteriors
(the purple contours) are pulled to lower $S_{8}$. As we have discussed
already, the direction is likely at least in part a function of the prior in
the $S_{8}-\Omega_{\mathrm{m}}$ plane, which restricts further movement in the
low $\Omega_{\mathrm{m}}$ direction. Interestingly the IA parameters are not
only tightened considerably but drawn towards the region favoured by $2\times
2$pt (the estimates of $A_{1}$, from cosmic shear + SR and $2\times 2$pt with
and without SR all agree very well, despite the information coming from very
different scales). The end result is something not dissimilar to the joint
$3\times 2$pt results shown in Section 4.1 \- the contours are tighter in the
$S_{8}$ direction and shifted, but are still not correctly centred.

The $2\times 2$pt contours on the other hand are stable, both in size and
position. Interestingly, the impact is even smaller than in Figure 6, where we
saw a $\sim 0.35\sigma$ shift when including scales in $\gamma_{t}$ down to
$4$ Mpc$/h$. We can draw a few things from this. First, as seen in previous
works (e.g. Sánchez, Prat et al. 2022; Amon et al. 2022) SR add most when IAs
(or redshift parameters) are not well constrained. In the $2\times 2$pt case,
the large scale data provide a reasonable constraint, and so there is less
impact in the $S_{8}-\Omega_{\mathrm{m}}$ plane. Second, we see that $A_{1}$
changes very little when the SR likelihood is included. This is despite the
fact that we are including scales where we expect the higher-order TATT
contributions to worsen. One possibility here is that the IA constraints are
simply dominated by the larger scales - this seems doubtful due to the
relatively high precision of the SR data, and the fact that adding SR pulls
the $1\times 2$pt constraints on $A_{1}$ towards the same value. Another is
that the TATT model on these scales can be consistently absorbed into the NLA
model, implying a relatively weak scale dependence across the range $2-6$
Mpc$/h$. Finally, even with the addition of SR we see only small $<1\sigma$
bias in $2\times 2$pt. We can take away from this that although we are
including more contaminated data vector, on scales where the S/N is relatively
high, the dynamic seen in a large scales only analysis is unchanged. That is,
the data can still distinguish fairly well between shifts in $S_{8}$ and
$\Omega_{\mathrm{m}}$ and unmodelled IA contributions. The inclusion of SR in
cosmic shear lessens the discrepancy slightly by shifting the cosmic shear
posterior, but we still see IA modelling error manifested as a tension between
probes.

We should finally note that this setup mirrors the choices made in DES Y3. One
could use SR purely as a redshift calibration tool by selecting bin pairs
where the ratios are geometric (see Sánchez, Prat et al. 2022 Figure 2 and
Section III). In this case, they will act only to tighten the posteriors on
the (already largely unbiased) redshift parameters. Given this, we would
expect even less impact compared with the solid contours in Figure 18 in such
an analysis. In short, our conclusions appear to be robust to the addition of
small scale shear ratios.


$\displaystyle|\langle\mathbf{E},\mathbf{Y}^{*}-\bar{\mathbf{Y}}\rangle|\leq\xi_{E}\left(\frac{\mu^{5}\kappa_{n}}{n}\right)^{1/2}\|\bar{\mathbf{Y}}-\mathbf{Y}^{*}\|_{F}.$
(59)

Now if $\|\bar{\mathbf{Y}}-\mathbf{Y}^{*}\|_{F}\geq
2\|\tilde{\mathbf{Y}}-\mathbf{Y}^{*}\|_{F}$, then
$\|\bar{\mathbf{Y}}-\mathbf{Y}^{*}\|_{F}\leq
2\|\bar{\mathbf{Y}}-\tilde{\mathbf{Y}}\|_{F}$. This, together with (58), (59),
and (52), yields

$\displaystyle\|\tilde{\mathbf{Y}}-\mathbf{Y}^{*}\|_{F}^{2}\leq$
$\displaystyle\|\tilde{\mathbf{Y}}-\bar{\mathbf{Y}}\|_{F}^{2}+4\xi_{E}\left(\frac{\mu^{5}\kappa_{n}}{n}\right)^{1/2}\|\bar{\mathbf{Y}}-\tilde{\mathbf{Y}}\|_{F}$
$\displaystyle\leq$ $\displaystyle
2\|\tilde{\mathbf{Y}}-\bar{\mathbf{Y}}\|_{F}^{2}+4\mu^{5}\frac{\kappa_{n}\xi_{E}^{2}}{n}$
$\displaystyle\leq$ $\displaystyle 16\mu
L_{0}^{2}\frac{nd_{n}\xi_{E}^{2}}{(n\delta_{n}-\xi_{E})^{2}}+4\mu^{5}\frac{\kappa_{n}\xi_{E}^{2}}{n}.$
(60)

On the other hand, if
$\|\bar{\mathbf{Y}}-\mathbf{Y}^{*}\|_{F}<2\|\tilde{\mathbf{Y}}-\mathbf{Y}^{*}\|_{F}$,
then (58) and (59) imply

$\displaystyle\|\tilde{\mathbf{Y}}-\mathbf{Y}^{*}\|_{F}^{2}\leq$
$\displaystyle\|\tilde{\mathbf{Y}}-\bar{\mathbf{Y}}\|_{F}^{2}+4\xi_{E}\left(\frac{\mu^{5}\kappa_{n}}{n}\right)^{1/2}\|\tilde{\mathbf{Y}}-\mathbf{Y}^{*}\|_{F}.$

By a similar argument used to obtain (55), this and (52) yield

$\displaystyle\|\tilde{\mathbf{Y}}-\mathbf{Y}^{*}\|_{F}^{2}\leq$
$\displaystyle
16\mu^{5}\frac{\kappa_{n}\xi_{E}^{2}}{n}+2\|\tilde{\mathbf{Y}}-\bar{\mathbf{Y}}\|_{F}^{2}$
$\displaystyle\leq$ $\displaystyle
16\mu^{5}\frac{\kappa_{n}\xi_{E}^{2}}{n}+16\mu
L_{0}^{2}\frac{nd_{n}\xi_{E}^{2}}{(n\delta_{n}-\xi_{E})^{2}}.$ (61)

In view of (56), (57), (60), (61) and (C5), the deisred result follows.

##  Appendix C Further technical details

In this section, we present some additional auxiliary results along with the
proofs of (36), (39), (40), (47), (52). Some existing results that are useful
in our proofs are also stated here for completeness with the references to
their proofs in the literature. These results are stated in the forms that are
most convenient for our use, which may not be in full generality.

######  Proposition 8 (Ruhe, 1970)

Let $\mathbf{A},\mathbf{B}$ be matrices with size $m\times n$ and $n\times p$
respectively. Then

$\displaystyle\sum_{j=1}^{n}\sigma_{j}^{2}(\mathbf{A})\sigma_{j}^{2}(\mathbf{B})\geq\|\mathbf{AB}\|_{F}^{2}\geq\sum_{j=1}^{n}\sigma_{n-j+1}^{2}(\mathbf{A})\sigma_{j}^{2}(\mathbf{B}).$

###### Remark 9

One consequence of this inequality we frequently use is
$\sigma_{1}^{2}(\mathbf{A})\|\mathbf{B}\|_{F}^{2}\geq\|\mathbf{AB}\|_{F}^{2}\geq\sigma_{n}^{2}(\mathbf{A})\|\mathbf{B}\|_{F}^{2}$.
Note also that by transposition the roles of $\mathbf{A}$ and $\mathbf{B}$ can
be interchanged on the left- and right-most expressions.

###### Lemma 10

Assume (C1)-(C2) and that $\sum_{j=1}^{p_{n}}\|\mathbf{B}_{j}^{*}\|_{*}\leq
L$. Suppose $L_{n}=d_{n}^{1/2}L_{0}$ is chosen so that $L_{0}\geq
L/(1-\epsilon_{L})$ with $1-\epsilon_{L}\leq 1/(4\mu^{2})$. Then for first-
and second-stage RGA, with probability tending to one,

(i)

$\displaystyle\inf_{k\geq
1}\frac{1}{nd_{n}}\|\mathbf{X}_{\hat{j}_{k}}\tilde{\mathbf{B}}_{\hat{j}_{k}}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}\geq$
$\displaystyle(1-\epsilon_{L})\mu L_{0}^{2}$ (62) $\displaystyle\inf_{k\geq
1}\frac{1}{nd_{n}}\|\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}\geq$
$\displaystyle(1-\epsilon_{L})\mu L_{0}^{2}$ (63)

(ii)

$\displaystyle\sup_{k\geq 1}|\lambda_{k}-\hat{\lambda}_{k}|\leq$
$\displaystyle\frac{2}{(1-\epsilon_{L})L_{0}}\frac{\xi_{E}}{n\sqrt{d_{n}}}$
(64)

(iii)

$\displaystyle\max_{1\leq k\leq K_{n}}\lambda_{k}\leq 1.$ (65)

Proof  We shall prove the results for the second-stage RGA. The corresponding
proofs for first-stage RGA follow similarly and thus are omitted. It is also
sufficient to prove (i)-(iii) assuming the condition described in (C1) holds
almost surely because the event that the condition holds has probability
tending to one. It will greatly simplify the exposition (without repeating
that the inequalities holds except on a vanishing event). Note that

$\displaystyle\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top},\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(k-1)}\rangle$
$\displaystyle=$
$\displaystyle\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top},\mathbf{Y}-\hat{\mathbf{G}}^{(k-1)}\rangle-\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top},\mathbf{E}\rangle$
$\displaystyle\geq$
$\displaystyle-|\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top},\mathbf{E}\rangle|$
$\displaystyle\geq$
$\displaystyle-\|\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{X}_{\hat{j}_{k}}^{\top}\mathbf{E}\|_{op}\|\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}\|_{*}$
$\displaystyle\geq$ $\displaystyle-\mu L_{n}\xi_{E},$

where the first inequality follows because
$\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{T},\mathbf{Y}-\hat{\mathbf{G}}^{(k-1)}\rangle\geq
0$ with probability one and the second inequality follows because the dual
norm of the nuclear norm is the operator norm. By Proposition 8, we have

$\displaystyle\|\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}$
$\displaystyle\geq$
$\displaystyle\|\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}\|_{F}^{2}-2\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top},\hat{\mathbf{G}}^{(k-1)}\rangle$
$\displaystyle\geq$ $\displaystyle
n\mu^{-1}\|\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}\|_{F}^{2}$
$\displaystyle+2\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top},\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(k-1)}\rangle-2\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top},\tilde{\mathbf{Y}}\rangle$
$\displaystyle\geq$ $\displaystyle n\mu^{-1}L_{n}^{2}-2\mu
L_{n}\xi_{E}-2\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top},\tilde{\mathbf{Y}}\rangle,$

where the last inequality follows from the fact that $\hat{\mathbf{S}}_{k}$ is
rank-one with singular value $L_{n}$. Thus, by writing
$\hat{\mathbf{S}}_{k}=L_{n}\mathbf{a}\mathbf{b}^{T}$ for some unit vectors
$\mathbf{a},\mathbf{b}$, we have
$\|\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{T}\|_{F}^{2}=L_{n}^{2}\|\mathbf{U}_{\hat{j}_{k}}\mathbf{a}\mathbf{b}^{T}\mathbf{V}_{\hat{j}_{k}}^{T}\|_{F}^{2}=L_{n}^{2}$.
Next, observe that

$\displaystyle|\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top},\tilde{\mathbf{Y}}\rangle|=$
$\displaystyle\left|\sum_{j=1}^{p_{n}}\langle\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top},\mathbf{X}_{j}\mathbf{B}_{j}^{*}\rangle\right|$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{p_{n}}\|\mathbf{B}_{j}^{*}\|_{*}\|\mathbf{X}_{j}^{\top}\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}\|_{op}$
$\displaystyle\leq$ $\displaystyle(1-\epsilon_{L})L_{n}^{2}n\mu.$

Therefore,

$\displaystyle(nd_{n})^{-1}\|\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}\geq$
$\displaystyle\mu^{-1}L_{0}^{2}-2(1-\epsilon_{L})L_{0}^{2}\mu-2\mu
L_{0}\frac{\xi_{E}}{n\sqrt{d_{n}}}$ $\displaystyle\geq$ $\displaystyle
2(1-\epsilon_{L})L_{0}^{2}\mu-2\mu L_{0}\frac{\xi_{E}}{n\sqrt{d_{n}}}.$

Since $\xi_{E}=o_{p}(n\sqrt{d_{n}})$ by (C2), (63) follows.

For (64), note first that if the solutions to the line search problems (9) and
(25) (with $\tilde{\mathbf{B}}_{\hat{j}_{k}}$ replaced by
$\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}$)
for second-stage RGA are not constrained to be in $[0,1]$, then they are given
by

$\displaystyle\hat{\lambda}_{k,uc}=$
$\displaystyle\frac{\langle\mathbf{Y}-\hat{\mathbf{G}}^{(k-1)},\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}-\hat{\mathbf{G}}^{(k-1)}\rangle}{\|\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}},$
$\displaystyle\lambda_{k,uc}=$
$\displaystyle\frac{\langle\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(k-1)},\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}-\hat{\mathbf{G}}^{(k-1)}\rangle}{\|\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}}.$

Since $\hat{\mathbf{G}}^{(l)}$ can always be expressed as
$\hat{\mathbf{G}}^{(l)}=\sum_{j\in\hat{J}}\mathbf{X}_{j}\hat{\mathbf{\Sigma}}_{j}^{-1}\mathbf{U}_{j}\mathbf{A}_{j}\mathbf{V}_{j}^{\top}$
with $\sum_{j\in\hat{J}}\|\mathbf{A}_{j}\|_{*}\leq L_{n}$, it follows that

$\displaystyle|\hat{\lambda}_{k}-\lambda_{k}|\leq|\hat{\lambda}_{k,uc}-\lambda_{k,uc}|=$
$\displaystyle\frac{|\langle\mathbf{E},\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}-\hat{\mathbf{G}}^{(k-1)}\rangle|}{\|\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}}$
$\displaystyle\leq$
$\displaystyle\frac{2L_{n}\mu\xi_{E}}{\|\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{\top}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}}$
$\displaystyle\leq$
$\displaystyle\frac{2\xi_{E}}{nd_{n}^{1/2}(1-\epsilon_{L})L_{0}},$

with probability tending to one, where the last inequality follows from (63).

For (65), it suffices to prove that
$\lim_{n\rightarrow\infty}\mathbb{P}(E_{n})=1$, where
$E_{n}=\left\\{\max_{1\leq k\leq K_{n}}\lambda_{k,uc}\leq 1\right\\}$. On
$E_{n}^{c}$, there exists some $k$ such that, by Cauchy-Schwarz inequality and
(27),

$\displaystyle\|\mathbf{X}_{\hat{j}_{k}}\hat{\mathbf{\Sigma}}_{\hat{j}_{k}}^{-1}\mathbf{U}_{\hat{j}_{k}}\hat{\mathbf{S}}_{k}\mathbf{V}_{\hat{j}_{k}}^{T}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}\leq$
$\displaystyle\|\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}$
$\displaystyle\leq$
$\displaystyle\|\tilde{\mathbf{Y}}\|_{F}^{2}+2\sum_{j=1}^{k-1}\langle\mathbf{E},\hat{\mathbf{G}}^{(k-j)}-\mathbf{G}^{(k-j)}\rangle$
$\displaystyle=$
$\displaystyle\|\tilde{\mathbf{Y}}\|_{F}^{2}+2\sum_{l=1}^{k-1}(\hat{\lambda}_{l}-\lambda_{l})\langle\mathbf{E},\mathbf{X}_{\hat{j}_{l}}\hat{\mathbf{\Sigma}}_{\hat{j}_{l}}^{-1}\mathbf{U}_{\hat{j}_{l}}\hat{\mathbf{S}}_{l}\mathbf{V}_{\hat{j}_{l}}^{\top}-\hat{\mathbf{G}}^{(l-1)}\rangle$
$\displaystyle\leq$
$\displaystyle\|\tilde{\mathbf{Y}}\|_{F}^{2}+4K_{n}L_{n}\mu\xi_{E}\max_{1\leq
l\leq k}|\hat{\lambda}_{l}-\lambda_{l}|.$ (66)

It is easy to see that

$\displaystyle\|\tilde{\mathbf{Y}}\|_{F}=\left\|\sum_{j=1}^{p_{n}}\mathbf{X}_{j}\mathbf{B}_{j}^{*}\right\|_{F}\leq(1-\epsilon_{L})L_{n}\sqrt{n\mu}.$
(67)

Thus, by (63), (64) and (C)-(67), we have

$\displaystyle\mathbb{P}(E_{n}^{c})\leq\mathbb{P}\left((1-\epsilon_{L})L_{0}^{2}\mu\\{1-(1-\epsilon_{L})\\}\leq\frac{8\mu}{1-\epsilon_{L}}\frac{K_{n}\xi_{E}^{2}}{n^{2}d_{n}}\right)+o(1)=o(1),$

where the last equality follows from (C2).

###### Lemma 11

Let $\\{a_{m}\\}$ be a nonnegative sequence of reals. If

$\displaystyle a_{0}\leq A,\mbox{ and }a_{m}\leq
a_{m-1}\left(1-\frac{\xi^{2}a_{m-1}}{A}\right)+b_{m},$

for $m=1,2,\ldots,$ where $b_{m}\geq 0$ with $b_{0}=0$, then for each $m$,

$\displaystyle a_{m}\leq\frac{A}{1+m\xi^{2}}+\sum_{k=0}^{m}b_{k}.$ (68)

Proof  We prove by induction. When $m=0$, (68) holds by assumption. Suppose
now that (68) holds for some $m\geq 1$. Then

$\displaystyle a_{m+1}\leq$ $\displaystyle
a_{m}\left(1-\frac{\xi^{2}a_{m}}{A}\right)+b_{m+1}$ $\displaystyle\leq$
$\displaystyle\frac{1}{a_{m}^{-1}+\xi^{2}/A}+b_{m+1}$ $\displaystyle\leq$
$\displaystyle\frac{1}{\left(\frac{A}{1+m\xi^{2}}+\sum_{k=0}^{m}b_{k}\right)^{-1}+\xi^{2}/A}+b_{m+1}$
$\displaystyle=$
$\displaystyle\frac{\frac{A}{1+m\xi^{2}}+\sum_{k=0}^{m}b_{k}}{1+\frac{\xi^{2}}{A}\left(\frac{A}{1+m\xi^{2}}+\sum_{k=0}^{m}b_{k}\right)}+b_{m+1}$
$\displaystyle\leq$
$\displaystyle\frac{A}{1+(m+1)\xi^{2}}+\sum_{k=0}^{m+1}b_{k},$

where the second inequality follows from $1-x\leq 1/(1+x)$ for $x\geq 0$.

###### Remark 12

Lemma 11 is a slight modification of Lemma 3.1 of Temlyakov (2000).

Proof [Proof of (36)] On $\mathcal{E}_{n}^{c}(m)$, there exists some $l\leq m$
such that

$\displaystyle\tilde{\tau}d_{n}^{1/2}\xi_{E}\geq\max_{\begin{subarray}{c}1\leq
j\leq p_{n}\\\ \|\mathbf{B}_{j}\|_{*}\leq
L_{n}\end{subarray}}\langle\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(l-1)},\mathbf{X}_{j}\mathbf{B}_{j}-\hat{\mathbf{G}}^{(l-1)}\rangle\geq\|\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(l-1)}\|_{F}^{2}.$

By (27) and Lemma 10(ii), it follows that, on $\mathcal{E}_{n}^{c}(m)$ except
for a vanishing event,

$\displaystyle\|\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(m)}\|_{F}^{2}\leq$
$\displaystyle\|\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(l-1)}\|_{F}^{2}+2\sum_{k=l}^{m}\langle\mathbf{E},\hat{\mathbf{G}}^{(k)}-\mathbf{G}^{(k)}\rangle$
$\displaystyle\leq$
$\displaystyle\tilde{\tau}d_{n}^{1/2}\xi_{E}+2\sum_{k=l}^{m}(\hat{\lambda}_{k}-\lambda_{k})\langle\mathbf{E},\mathbf{X}_{\hat{j}_{k}}\tilde{\mathbf{B}}_{\hat{j}_{k}}-\hat{\mathbf{G}}^{(k-1)}\rangle$
$\displaystyle\leq$
$\displaystyle\tilde{\tau}d_{n}^{1/2}\xi_{E}+\frac{8m\xi_{E}^{2}}{n(1-\epsilon_{L})},$

which is the desired result.

Proof [Proof of (39) and (40)] Note first that for any $D>0$, $(D+x)/(D-x)\leq
1+3x/D$ for all $0\leq x\leq(1-\sqrt{2/3})D$. It is not difficult to see that

$\displaystyle\mathbb{P}\left\\{\frac{4L_{0}\xi_{E}}{nd_{n}^{1/2}}\leq(1-\sqrt{\frac{2}{3}})\left((nd_{n})^{-1}\|\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(k)}\|_{F}^{2}+(nd_{n})^{-1}\|\mathbf{E}\|_{F}^{2}\right),1\leq
k\leq\hat{k},\mathcal{G}_{n}\right\\}$ $\displaystyle\geq$
$\displaystyle\mathbb{P}\left\\{\frac{4L_{0}\xi_{E}}{nd_{n}^{1/2}}\leq(1-\sqrt{\frac{2}{3}})M^{-1}\right\\}-o(1)$
$\displaystyle\rightarrow$ $\displaystyle 1.$

Thus, on $\mathcal{G}_{n}$ except for a vanishing event,

$\displaystyle A_{k}\leq$ $\displaystyle
1+\frac{12L_{0}\xi_{E}/(nd_{n}^{1/2})}{(nd_{n})^{-1}\|\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(k)}\|_{F}^{2}+(nd_{n})^{-1}\|\mathbf{E}\|_{F}^{2}}$
$\displaystyle\leq$ $\displaystyle 1+12ML_{0}\frac{\xi_{E}}{nd_{n}^{1/2}},$

for all $1\leq k\leq\hat{k}$. This proves (39). We now turn to (40). Since for
any positive $A$ and $B$, $A/(B+x)\geq A(1-x/B)/B$ for all $x\geq 0$, it
follows from (37) that on $\mathcal{G}_{n}$ except for a vanishing event,

$\displaystyle B_{k}\geq$
$\displaystyle\frac{\tau^{2}s_{n}^{-1}}{4L_{0}^{2}\mu^{2}}\frac{(nd_{n})^{-1}\|\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}}{(nd_{n})^{-1}\|\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}+(nd_{n})^{-1}\|\mathbf{E}\|_{F}^{2}}$
$\displaystyle\times\left(1-\frac{4L_{0}\xi_{E}/(nd_{n}^{1/2})}{(nd_{n})^{-1}\|\tilde{\mathbf{Y}}-\hat{\mathbf{G}}^{(k-1)}\|_{F}^{2}+(nd_{n})^{-1}\|\mathbf{E}\|_{F}^{2}}\right)$
$\displaystyle\geq$
$\displaystyle\frac{\tau^{2}s_{n}^{-1}}{4L_{0}^{2}\mu^{2}}\frac{1}{1+\mu
Ms_{n}}\left(1-\frac{4ML_{0}\xi_{E}}{nd_{n}^{1/2}}\right)$

for $1\leq k\leq\hat{k}$, which proves (40).

Proof [Proof of (47)] Let

$\displaystyle\mathbf{H}=\sum_{j\in\hat{J}}\mathbf{X}_{j}\hat{\mathbf{\Sigma}}_{j}^{-1}\mathbf{U}_{j}\mathbf{D}_{j}\mathbf{V}_{j}^{\top}\in\mathcal{B}.$

Note that Proposition 8 and (C3) imply

$\displaystyle\|\bar{\mathbf{Y}}-\mathbf{H}\|_{F}^{2}\geq$ $\displaystyle
n\mu^{-1}\left\\{\sum_{j\in\hat{J}_{o}}\|\hat{\mathbf{\Sigma}}_{j}^{-1}\mathbf{U}_{j}(\mathbf{L}_{j}\mathbf{\Lambda}_{j}\mathbf{R}_{j}^{\top}-\mathbf{D}_{j})\mathbf{V}_{j}^{\top}\|_{F}^{2}+\sum_{j\in\hat{J}-\hat{J}_{o}}\|\hat{\mathbf{\Sigma}}_{j}^{-1}\mathbf{U}_{j}\mathbf{D}_{j}\mathbf{V}_{j}^{\top}-\mathbf{B}_{j}^{*}\|_{F}^{2}\right\\}$
$\displaystyle\geq$ $\displaystyle
n\mu^{-3}\left\\{\sum_{j\in\hat{J}_{o}}\|\mathbf{L}_{j}\mathbf{\Lambda}_{j}\mathbf{R}_{j}^{\top}-\mathbf{D}_{j}\|_{F}^{2}+\sum_{j\in\hat{J}-\hat{J}_{o}}\|\mathbf{U}_{j}^{\top}\hat{\mathbf{\Sigma}}_{j}\mathbf{B}_{j}^{*}\mathbf{V}_{j}-\mathbf{D}_{j}\|_{F}^{2}\right\\}$
$\displaystyle\geq$
$\displaystyle\frac{n}{\mu^{3}\kappa_{n}}\left\\{\sum_{j\in\hat{J}_{o}}\|\mathbf{L}_{j}\mathbf{\Lambda}_{j}\mathbf{R}_{j}^{\top}-\mathbf{D}_{j}\|_{*}+\sum_{j\in\hat{J}-\hat{J}_{o}}\|\mathbf{U}_{j}^{\top}\hat{\mathbf{\Sigma}}_{j}\mathbf{B}_{j}^{*}\mathbf{V}_{j}-\mathbf{D}_{j}\|_{*}\right\\}^{2}.$

Since $\mathbf{H}\in\mathcal{B}$, we have

$\displaystyle\left\\{\sum_{j\in\hat{J}_{o}}\|\mathbf{L}_{j}\mathbf{\Lambda}_{j}\mathbf{R}_{j}^{\top}-\mathbf{D}_{j}\|_{*}+\sum_{j\in\hat{J}-\hat{J}_{o}}\|\mathbf{U}_{j}^{\top}\hat{\mathbf{\Sigma}}_{j}\mathbf{B}_{j}^{*}\mathbf{V}_{j}-\mathbf{D}_{j}\|_{*}\right\\}^{2}\leq\frac{9d_{n}L_{0}^{2}}{16}=\frac{9L_{n}^{2}}{16}.$

By the triangle inequality, we have
$\sum_{j\in\hat{J}}\|\mathbf{D}_{j}\|_{*}\leq
3L_{n}/4+\sum_{j\in\hat{J}_{o}}\|\mathbf{\Lambda}_{j}\|_{*}+\sum_{j\in\hat{J}-\hat{J}_{o}}\|\hat{\mathbf{\Sigma}}_{j}\mathbf{B}_{j}^{*}\|_{*}$.
Because of (C6), and $\hat{J}_{o}\subset J_{o}$ (with probability tending to
one),
$\sum_{j\in\hat{J}_{o}}\|\mathbf{\Lambda}_{j}\|_{*}+\sum_{j\in\hat{J}-\hat{J}_{o}}\|\hat{\mathbf{\Sigma}}_{j}\mathbf{B}_{j}^{*}\|_{*}\leq\sum_{j\in\hat{J}}\|\hat{\mathbf{\Sigma}}_{j}\mathbf{B}_{j}^{*}\|_{*}\leq\mu(1-\epsilon_{L})L_{n}\leq
4^{-1}\mu^{-1}L_{n}\leq L_{n}/4$. Hence
$\sum_{j\in\hat{J}_{\hat{k}}}\|\mathbf{D}_{j}\|_{*}\leq L_{n}$, which proves
$\mathbf{H}\in\mathcal{C}_{L}$.

###### Proposition 13

Let $\mathbf{A}^{*}$ be an $m\times n$ matrix and
$\mathbf{A}=\mathbf{A}^{*}+\mathbf{E}$ be its perturbed version. Let
$\mathbf{U}_{*}\mathbf{\Sigma}_{*}\mathbf{V}_{*}^{\top}$ and
$\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ be their truncated SVD of rank
$r_{*}$, respectively. If
$\sigma_{r_{*}}(\mathbf{A}^{*}):=\sigma_{r_{*}}>\sigma_{r_{*}+1}(\mathbf{A}^{*})=0$,
and if $\|\mathbf{E}\|_{op}<\sigma_{r_{*}}$, then

$\displaystyle\max\\{\mathrm{dist}(\mathbf{U}_{*},\mathbf{U}),\mathrm{dist}(\mathbf{V}_{*},\mathbf{V})\\}\leq\frac{\sqrt{2}\max\\{\|\mathbf{E}^{\top}\mathbf{U}_{*}\|_{op},\|\mathbf{E}\mathbf{V}_{*}\|_{op}\\}}{\sigma_{r_{*}}-\|\mathbf{E}\|_{op}},$

where
$\mathrm{dist}(\mathbf{Q},\mathbf{Q}_{*})=\min_{\mathbf{R}}\|\mathbf{QR}-\mathbf{Q}_{*}\|_{op}$
for any two orthogonal matrices $\mathbf{Q}$, $\mathbf{Q}^{*}$ with $r$
columns, where the minimum is taken over all $r\times r$ orthonormal matrices.

###### Remark 14

Proposition 13 is a consequence of the perturbation bounds for singular values
(Wedin, 1972). A proof can be found in Chen et al. (2021).

Proof [Proof of (52)] Note first that

$\displaystyle\bar{\mathbf{Y}}-\tilde{\mathbf{Y}}=$
$\displaystyle\sum_{j\in\hat{J}_{o}}\mathbf{X}_{j}\hat{\mathbf{\Sigma}}_{j}^{-1}(\mathbf{U}_{j}\mathbf{L}_{j}-\tilde{\mathbf{U}}_{j})\mathbf{\Lambda}_{j}\tilde{\mathbf{V}}_{j}^{\top}$
$\displaystyle+\sum_{j\in\hat{J}_{o}}\mathbf{X}_{j}\hat{\mathbf{\Sigma}}_{j}^{-1}\mathbf{U}_{j}\mathbf{L}_{j}\mathbf{\Lambda}_{j}(\mathbf{V}_{j}\mathbf{R}_{j}-\tilde{\mathbf{V}}_{j})^{\top}.$

By triangle inequality,

$\displaystyle\|\bar{\mathbf{Y}}-\tilde{\mathbf{Y}}\|_{F}\leq$
$\displaystyle\sqrt{n\mu}\left(\sum_{j\in\hat{J}_{o}}\|\mathbf{\Lambda}_{j}\|_{F}\right)\left\\{\max_{j\in\hat{J}_{o}}\|\mathbf{U}_{j}\mathbf{L}_{j}-\tilde{\mathbf{U}}_{j}\|_{op}+\max_{j\in\hat{J}_{o}}\|\mathbf{V}_{j}\mathbf{R}_{j}-\tilde{\mathbf{V}}_{j}\|_{op}\right\\}.$
(69)

Let $\mathbf{U}_{j,\bar{r}_{j}}$ and $\mathbf{V}_{j,\bar{r}_{j}}$ be sub-
matrices of $\mathbf{U}_{j}$ and $\mathbf{V}_{j}$ consisting of column vectors
that correspond to the leading $\bar{r}_{j}$ singular vectors. Write
$\mathbf{U}_{j}=(\mathbf{U}_{j,\bar{r}_{j}},\mathbf{U}_{j,-\bar{r}_{j}})$ and
$\mathbf{V}_{j}=(\mathbf{V}_{j,\bar{r}_{j}},\mathbf{V}_{j,-\bar{r}_{j}})$.
Since
$\mathbf{X}_{j}^{\top}\tilde{\mathbf{Y}}=\mathbf{X}_{j}^{\top}\mathbf{Y}-\mathbf{X}_{j}^{\top}\mathbf{E}$,
it follows from Proposition 13 and (C5) that there exist
$\bar{r}_{j}\times\bar{r}_{j}$ orthonormal matrices $\tilde{\mathbf{L}}_{j}$
and $\tilde{\mathbf{R}}_{j}$ such that with probability tending to one,

$\displaystyle\max\left\\{\|\mathbf{U}_{j,\bar{r}_{j}}\tilde{\mathbf{L}}_{j}-\tilde{\mathbf{U}}_{j}\|_{op},\|\mathbf{V}_{j,\bar{r}_{j}}\tilde{\mathbf{R}}_{j}-\tilde{\mathbf{V}}_{j}\|_{op}\right\\}\leq$
$\displaystyle\frac{\sqrt{2}\max\\{\|\mathbf{E}^{\top}\mathbf{X}_{j}\tilde{\mathbf{U}}_{j}\|_{op},\|\mathbf{X}_{j}^{\top}\mathbf{E}\tilde{\mathbf{V}}_{j}\|_{op}\\}}{n\delta_{n}-\|\mathbf{X}_{j}^{\top}\mathbf{E}\|_{op}}$
$\displaystyle\leq$
$\displaystyle\frac{\sqrt{2}\xi_{E}}{n\delta_{n}-\xi_{E}}.$

Set
$\mathbf{L}_{j}^{\top}=(\tilde{\mathbf{L}}_{j}^{\top},\mathbf{0}_{\bar{r}_{j}\times(\hat{r}-\bar{r}_{j})})$
and
$\mathbf{R}_{j}^{\top}=(\tilde{\mathbf{R}}_{j}^{\top},\mathbf{0}_{\bar{r}_{j}\times(\hat{r}-\bar{r}_{j})})$
for $j\in\hat{J}_{o}$ in (69). Then by (C4) and (C6), it follows that

$\displaystyle\|\bar{\mathbf{Y}}-\tilde{\mathbf{Y}}\|_{F}^{2}\leq
n\mu\left(\sum_{j\in\hat{J}_{o}}\|\mathbf{\Lambda}_{j}\|_{F}\right)^{2}\left(\frac{2\sqrt{2}\xi_{E}}{n\delta_{n}-\xi_{E}}\right)^{2}\leq
8\mu L^{2}nd_{n}\frac{\xi_{E}^{2}}{(n\delta_{n}-\xi_{E})^{2}}.$

Proof [Proof of Corollary 5] By Lemma 2,
$\sharp(\hat{J})+\hat{r}=O_{p}(s_{n}^{2})$. Thus running the first-stage RGA
with the just-in-time stopping criterion costs

$\displaystyle O_{p}(s_{n}^{2}(n_{1}+d_{n}))$ (70)

bytes of communication per computing node. In addition, preparing
$\\{\hat{\mathbf{\Sigma}}_{j}^{-1}:j\in\hat{J}\\}$ and
$(\mathbf{U}_{j},\mathbf{V}_{j})$ for $j\in\hat{J}$ with $q_{n,j}\wedge
d_{n}>\hat{r}$ costs

$\displaystyle
O_{p}\left(\sum_{j\in\hat{J}}\\{q_{n,j}^{2}+(q_{n,j}d_{n}+\hat{r}(q_{n,j}+d_{n}))\mathbf{1}\\{q_{n,j}\wedge
d_{n}>\hat{r}\\}\\}\right)$ $\displaystyle=$ $\displaystyle
O_{p}(n_{1}^{2\alpha}s_{n}^{2}+n_{1}^{\alpha}d_{n}s_{n}^{2}+s_{n}^{4}(n_{1}^{\alpha}+d_{n})).$
(71)

Since the communication costs per node at the $k$-th iteration of the second-
stage RGA is at most

$\displaystyle
O_{p}\left(\sum_{j\in\hat{J}}\left(\hat{r}^{2}\mathbf{1}\\{q_{n,j}\wedge
d_{n}>\hat{r}\\}+q_{n,j}d_{n}\mathbf{1}\\{q_{n,j}\wedge
d_{n}\leq\hat{r}\\}\right)+d_{n}k+n_{1}\right)$ $\displaystyle=$
$\displaystyle
O_{p}\left(s_{n}^{6}+n_{1}^{\alpha}d_{n}s_{n}^{2}+d_{n}k+n_{1}\right),$

running $m_{n}=O_{p}(s_{n}^{4}\log(n^{2}d_{n}/\xi_{n}^{2}))$ iterations (see
Theorem 3 for the definition of $m_{n}$) costs

$\displaystyle
O_{p}\left((s_{n}^{6}+s_{n}^{2}n_{1}^{\alpha}d_{n}+n_{1})s_{n}^{4}\log\frac{n^{2}d_{n}}{\xi_{n}^{2}}+d_{n}s_{n}^{8}\left(\log\frac{n^{2}d_{n}}{\xi_{n}^{2}}\right)^{2}\right).$
(72)

Combining (70)-(72) yields the desired result.

##  Appendix D TSRGA for high-dimensional generalized linear models

In this section, we apply the idea of TSRGA to and propose a modified
algorithm for estimating the generalized linear model (GLM). Focusing on the
case of a scalar response $y_{t}$, the GLM postulates that the probability
density function $f$ of $y_{t}$ (or the probability mass function if $y_{t}$
is discrete) belongs to the exponential family. In particular,

$\displaystyle f(y;\theta)=\exp[y\theta-r(\theta)+h(y)],$

and

$\displaystyle\mathbb{E}(y_{t}|{x}_{t,1},\ldots,x_{t,p_{n}})=r^{\prime}\left(\sum_{j=1}^{p_{n}}\beta_{j}^{*}x_{t,j}\right)$

where $\theta$ is called the natural parameter; $r$, $h$ are known functions,
and $r^{\prime}$ is the derivative of $r$, which is also known as the inverse
of the link function (see, e.g., Dunn and Smyth, 2018; Han et al., 2023). To
maximize the log-likelihood function, scaled as $y\theta-r(\theta)$, one can
minimize the following loss function

$\displaystyle\mathcal{L}_{n}(\mathbf{X}\bm{\beta})=\frac{1}{n}\sum_{t=1}^{n}\left[-y_{t}\left(\sum_{j=1}^{p_{n}}\beta_{j}x_{t,j}\right)+r\left(\sum_{j=1}^{p_{n}}\beta_{j}x_{t,j}\right)\right],$

where
$\mathcal{L}_{n}(\bm{\tau})=n^{-1}\sum_{t=1}^{n}(y_{t}\tau_{t}-r(\tau_{t}))$
for $\bm{\tau}=(\tau_{1},\ldots,\tau_{n})^{\top}$.

Interpreting $y_{t}-r^{\prime}(\sum_{j=1}^{p_{n}}\beta_{j}x_{t,j})$ as the
residual, we can implement RGA as follows. First initialize
$\hat{\mathbf{G}}^{(0)}=0$. Then for $k=1,2,\ldots,K_{n}$, find

$\displaystyle\hat{j}_{k}\in\arg\max_{1\leq j\leq
p_{n}}\left|\frac{1}{n}\sum_{t=1}^{n}\left(y_{t}-r^{\prime}(\hat{G}^{(k-1)}_{t})\right){x}_{t,j}\right|$
(73)

and update

$\displaystyle\hat{\mathbf{G}}^{(k)}=(1-\hat{\lambda}_{k})\hat{\mathbf{G}}^{(k-1)}+\hat{\lambda}_{k}Ls_{k}\mathbf{z}_{\hat{j}_{k}},$
(74)

where
$\hat{\mathbf{G}}^{(k)}=(\hat{G}_{1}^{(k)},\ldots,\hat{G}_{n}^{(k)})^{\top}$,
$L>0$ is given, $\mathbf{z}_{j}=(x_{1,j},\ldots,x_{n,j})^{\top}$,

$\displaystyle
s_{k}=\mathrm{sgn}\left(\frac{1}{n}\sum_{t=1}^{n}\left(y_{t}-r^{\prime}(\hat{G}^{(k-1)}_{t})\right){x}_{t,\hat{j}_{k}}\right),$

and $\hat{\lambda}_{k}$ is determined by

$\displaystyle\hat{\lambda}_{k}=\arg\min_{\lambda\in[0,1]}\mathcal{L}_{n}((1-\lambda)\hat{\mathbf{G}}^{(k-1)}+\lambda
Ls_{k}\mathbf{z}_{\hat{j}_{k}}).$

It is not difficult to see that (73) can be easily solved for feature-
distributed data and constructing $\hat{\mathbf{G}}^{(k)}$ in each node
requires a communication cost of $O(n)$ bytes. The second-stage RGA can be
implemented similarly with the set of predictors considered in (73) restricted
to $\hat{J}$, the set of predictors chosen by the first-stage when the just-
in-time criterion is met. Finally, since $\mathcal{L}_{n}$ could take negative
values, we modify the just-in-time criterion (6) as

$\displaystyle\hat{k}=\min\left\\{1\leq k\leq
K_{n}:\left|\frac{\mathcal{L}_{n}(\hat{\mathbf{G}}^{(k)})}{\mathcal{L}_{n}(\hat{\mathbf{G}}^{(k-1)})}-1\right|<t_{n}\right\\}.$
(75)

In the same spirit as (6), (75) terminates the first-stage RGA as soon as the
improvement in the loss function is below certain threshold, which would save
some communication costs and speed up the algorithm.

Next, we examine the performance of this version of TSRGA ((73)-(75)) using
simulations. In the following experiments, the predictors $x_{t,j}$ are
generated as in Specification 2. We consider the following two specifications.

###### Specification 5

(Logit model) The response $y_{t}$ takes only values in $\\{0,1\\}$ and is
generated via

$\displaystyle\mathbb{P}(y_{t}=y,\theta_{t})=\theta_{t}^{y}(1-\theta_{t})^{1-y},\quad\theta_{t}=\frac{1}{1+\exp(-\sum_{j=1}^{p_{n}}\beta_{j}^{*}x_{t,j})}$

where
$(\beta_{1}^{*},\beta_{2}^{*},\beta_{3}^{*},\beta_{4}^{*},\beta_{5}^{*})=(-2.4,1.8,-1.9,2.8,-2.2)$,
$\beta_{j}^{*}=0$ for $j>5$. For this model, we have
$r(\theta)=\log(1+\exp(\theta))$.

###### Specification 6

(Poisson model) The response $y_{t}$ takes values in $\\{0,1,2,\ldots\\}$ and
is generated via

$\displaystyle\mathbb{P}(y_{t}=y,\theta_{t})=\frac{\theta_{t}^{y}e^{-\theta_{t}}}{y!},\quad\theta_{t}=\exp(\sum_{j=1}^{p_{n}}\beta_{j}^{*}x_{t,j})$

and
$(\beta_{1}^{*},\beta_{2}^{*},\beta_{3}^{*},\beta_{4}^{*},\beta_{5}^{*})=(0.15,-0.25,0.35,-0.45,0.55)$,
$\beta_{j}^{*}=0$ for $j>5$. For this model, we have $r(\theta)=\exp(\theta)$.

As a benchmark, we compare with the $\ell_{1}$-regularized GLM which solves

$\displaystyle\min_{\bm{\beta}}\mathcal{L}_{n}(\mathbf{X}\bm{\beta})+\lambda\|\bm{\beta}\|_{1}$
(76)

with $\lambda$ selected by 5-fold cross validation. Table 4 reports the
parameter estimation error $\|\hat{\bm{\beta}}-\bm{\beta}^{*}\|_{2}$, the
number of irrelevant variables selected (false positives, FP) and the number
of relevant variables not selected (false negatives, FN). For the logit model,
we additionally report the out-of-sample prediction accuracy on a test set of
size 500. For the Poisson model, the out-of-sample prediction error is
measured by RMSE. All these figures are averages over 500 independent
simulations.

Logit | $n=800,p=1200$ | $n=1200,p=2000$
---|---|---
| TSRGA |  $\ell_{1}$-GLM | TSRGA |  $\ell_{1}$-GLM
$\|\hat{\bm{\beta}}-\bm{\beta}^{*}\|_{2}$ | 0.698 | 2.185 | 0.689 | 2.036
FP | 0.018 | 82.808 | 0 | 105.070
FN | 0 | 0 | 0 | 0
Accuracy | 0.901 | 0.888 | 0.901 | 0.892
Poisson |  |  |  |
$\|\hat{\bm{\beta}}-\bm{\beta}^{*}\|_{2}$ | 0.135 | 0.190 | 0.060 | 0.144
FP | 6.638 | 25.470 | 1.830 | 25.146
FN | 0.114 | 0.008 | 0.020 | 0
RMSE | 1.324 | 1.363 | 1.283 | 1.329
Table 4: Simulation results for estimating high-dimensional GLMs.
$\ell_{1}$-GLM is defined in (76). The results are based on 500 simulations.

The results show that for both the logit and Poisson models, TSRGA yields
parsimonious and accurate coefficient estimates, with comparable out-of-sample
prediction accuracy to the $\ell_{1}$-GLM defined by (76). In particular, the
low FP and FN values of TSRGA may be due to its variable selection properties.
Though we expect the general conclusions about TSRGA in this paper, such as
the sure-screening property, to hold under the GLM framework, the rigorous
mathematical treatment is left for future work.

##  Appendix E Complementary simulation results

In this section, we present some additional simulation results regarding
Specifications 1 and 2. Figures 5 and 6 plot the parameter estimation error,
as in Figures 1 and 2, against the elapsed time. Clearly, TSRGA converges
within the least amount of time. In particular, its second-stage only takes a
very short amount of time, thanks to the dimension reduction after the just-
in-time stopping criterion. Other methods behave similarly as those in Figures
1 and 2, as their implementation cost scales directly with the number of
iterations.

Figures 7 and 8 plot the out-of-sample prediction error (measures by the root
mean square prediction error on an independent test sample) of the methods
under Specifications 1 and 2. For Specification 1, the final prediction
accuracy of TSRGA, cross-validated Lasso, and Hydra are similar. However, for
Specification 2, TSRGA clearly is the most desirable prediction tool among the
methods under consideration.

(a) $n=800$, $p_{n}=1200$

(b) $n=1200$, $p_{n}=2000$

(c) $n=1500$, $p_{n}=3000$

Figure 5: Logarithm of parameter estimation errors of various methods against
the elapsed time under Specification 1, where $n$ is the sample size and
$p_{n}$ is the dimension of predictors. The results are based on 100
simulations.

(a) $n=800$, $p_{n}=1200$

(b) $n=1200$, $p_{n}=2000$

(c) $n=1500$, $p_{n}=3000$

Figure 6: Logarithm of parameter estimation errors of various methods against
the elapsed time under Specification 2, where $n$ is the sample size and
$p_{n}$ is the dimension of predictors. The results are based on 100
simulations.

(a) $n=800$, $p_{n}=1200$

(b) $n=1200$, $p_{n}=2000$

(c) $n=1500$, $p_{n}=3000$

Figure 7: Logarithm of out-of-sample prediction errors of various methods
under Specification 1, where $n$ is the sample size and $p_{n}$ is the
dimension of predictors. The results are based on 100 simulations.

(a) $n=800$, $p_{n}=1200$

(b) $n=1200$, $p_{n}=2000$

(c) $n=1500$, $p_{n}=3000$

Figure 8: Logarithm of out-of-sample prediction errors of various methods
under Specification 2, where $n$ is the sample size and $p_{n}$ is the
dimension of predictors. The results are based on 100 simulations.

## References

  * Bellet et al. (2015) Aurélien Bellet, Yingyu Liang, Alireza Bagheri Garakani, Maria-Florina Balcan, and Fei Sha.  A distributed frank-wolfe algorithm for communication-efficient sparse learning.  In _Proceedings of the 2015 SIAM International Conference on Data Mining_ , pages 478–486, 2015.
  * Bertrand and Moonen (2010) A. Bertrand and M. Moonen.  Distributed adaptive node-specific signal estimation in fully connected sensor networks—part i: Sequential node updating.  _IEEE Transactions on Signal Processing_ , 58(10):5277–5291, 2010.
  * Bertrand and Moonen (2015) A. Bertrand and M. Moonen.  Distributed canonical correlation analysis in wireless sensor networks with application to distributed blind source separation.  _IEEE Transactions on Signal Processing_ , 63(18):4800–4813, 2015.
  * Bertrand and Moonen (2014) Alexander Bertrand and Marc Moonen.  Distributed adaptive estimation of covariance matrix eigenvectors in wireless sensor networks with application to distributed pca.  _Signal Processing_ , 104:120–135, 2014.
  * Bunea et al. (2011) Florentina Bunea, Yiyuan She, and Marten H. Wegkamp.  Optimal selection of reduced rank estimators of high-dimensional matrices.  _The Annals of Statistics_ , 39(2):1282–1309, 2011.
  * Bybee et al. (2021) Leland Bybee, Bryan T Kelly, Asaf Manela, and Dacheng Xiu.  Business news and business cycles.  Working Paper 29344, National Bureau of Economic Research, October 2021\.
  * Chen et al. (2013) Kun Chen, Hongbo Dong, and Kung-Sik Chan.  Reduced rank regression via adaptive nuclear norm penalization.  _Biometrika_ , 100(4):901–920, 09 2013.
  * Chen et al. (2021) Yuxin Chen, Yuejie Chi, Jianqing Fan, and Cong Ma.  Spectral methods for data science: A statistical perspective.  _Foundations and Trends®in Machine Learning_ , 14(5):566–806, 2021.
  * Dalcín and Fang (2021) Lisandro Dalcín and Yao-Lung L. Fang.  mpi4py: Status update after 12 years of development.  _Computing in Science and Engineering_ , 23(4):47–54, 2021.
  * Dalcín et al. (2005) Lisandro Dalcín, Rodrigo Paz, and Mario Storti.  Mpi for python.  _Journal of Parallel and Distributed Computing_ , 65(9):1108–1115, 2005.
  * Deerwester et al. (1990) Scott Deerwester, Susan T. Dumais, George W. Furnas, Thomas K. Landauer, and Richard Harshman.  Indexing by latent semantic analysis.  _Journal of the American Society for Information Science_ , 41(6):391–407, 1990.
  * Ding et al. (2020) Lijun Ding, Yingjie Fei, Qiantong Xu, and Chengrun Yang.  Spectral frank-wolfe algorithm: Strict complementarity and linear convergence.  In _Proceedings of the 37th International Conference on Machine Learning_ , volume 119, pages 2535–2544, 13–18 Jul 2020.
  * Ding et al. (2021) Lijun Ding, Jicong Fan, and Madeleine Udell.  $k$ fw: A frank-wolfe style algorithm with stronger subproblem oracles.  _arXiv preprint arXiv:2006.16142_ , 2021.
  * Dunn and Smyth (2018) P.K. Dunn and G.K. Smyth.  _Generalized Linear Models With Examples in R_.  Springer Texts in Statistics. Springer New York, 2018.
  * Fan and Lv (2008) Jianqing Fan and Jinchi Lv.  Sure independence screening for ultrahigh dimensional feature space.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 70(5):849–911, 2008.
  * Fan et al. (2015) Yingying Fan, Gareth M. James, and Peter Radchenko.  Functional additive regression.  _The Annals of Statistics_ , 43(5):2296–2325, 2015.
  * Fercoq et al. (2014) Olivier Fercoq, Zheng Qu, Peter Richtárik, and Martin Takáč.  Fast distributed coordinate descent for non-strongly convex losses.  In _2014 IEEE International Workshop on Machine Learning for Signal Processing (MLSP)_ , pages 1–6, 2014.
  * Frank and Wolfe (1956) Marguerite Frank and Philip Wolfe.  An algorithm for quadratic programming.  _Naval Research Logistics Quarterly_ , 3:95–110, 1956.
  * Gandhi et al. (2019) Priyank Gandhi, Tim Loughran, and Bill McDonald.  Using annual report sentiment as a proxy for financial distress in u.s. banks.  _Journal of Behavioral Finance_ , 20(4):424–436, 2019.
  * Gao and Tsay (2022) Zhaoxing Gao and Ruey S. Tsay.  Divide-and-conquer: A distributed hierarchical factor approach to modeling large-scale time series data.  _Journal of the American Statistical Association_ , 2022.
  * Garber (2020) Dan Garber.  Revisiting frank-wolfe for polytopes: Strict complementarity and sparsity.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 18883–18893, 2020.
  * Gene H. Golub and Wahba (1979) Michael Heath Gene H. Golub and Grace Wahba.  Generalized cross-validation as a method for choosing a good ridge parameter.  _Technometrics_ , 21(2):215–223, 1979.
  * Han et al. (2023) Yuefeng Han, Ruey S. Tsay, and Wei Biao Wu.  High dimensional generalized linear models for temporal dependent data.  _Bernoulli_ , 29(1):105–131, 2023.
  * Hanley and Hoberg (2019) Kathleen Weiss Hanley and Gerard Hoberg.  Dynamic interpretation of emerging risks in the financial sector.  _The Review of Financial Studies_ , 32(12):4543–4603, 02 2019.
  * Heinze et al. (2016) Christina Heinze, Brian McWilliams, and Nicolai Meinshausen.  Dual-loco: Distributing statistical estimation using random projections.  In _Proceedings of the 19th International Conference on Artificial Intelligence and Statistics_ , volume 51, pages 875–883, Cadiz, Spain, 2016.
  * Hu et al. (2019) Yaochen Hu, Di Niu, Jianming Yang, and Shengping Zhou.  Fdml: A collaborative machine learning framework for distributed features.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 2232–2240, 2019.
  * Ing (2020) Ching-Kang Ing.  Model selection for high-dimensional linear regression with dependent observations.  _The Annals of Statistics_ , 48(4):1959–1980, 2020.
  * Ing and Lai (2011) Ching-Kang Ing and Tze Leung Lai.  A stepwise regression method and consistent model selection for high-dimensional sparse linear models.  _Statistica Sinica_ , 21(4):1473–1513, 2011.
  * Jaggi (2013) Martin Jaggi.  Revisiting frank-wolfe: Projection-free sparse convex optimization.  _Proceedings of the 30th International Conference on Machine Learning_ , 28(1):427–435, 2013.
  * Jaggi and Lacoste-Julien (2015) Martin Jaggi and Simon Lacoste-Julien.  On the global linear convergence of frank-wolfe optimization variants.  _Advances in Neural Information Processing Systems_ , 28, 2015.
  * Jegadeesh and Wu (2013) Narasimhan Jegadeesh and Di Wu.  Word power: A new approach for content analysis.  _Journal of Financial Economics_ , 110(3):712–729, 2013.
  * Kogan et al. (2009) Shimon Kogan, Dimitry Levin, Bryan R Routledge, Jacob S Sagi, and Noah A Smith.  Predicting risk from financial reports with regression.  In _Proceedings of Human Language Technologies: The 2009 Annual Conference of the North American Chapter of the Association for Computational Linguistics_ , pages 272–280, 2009.
  * Lei et al. (2019) Qi Lei, Jiacheng Zhuo, Constantine Caramanis, Inderjit S Dhillon, and Alexandros G Dimakis.  Primal-dual block generalized frank-wolfe.  In _Advances in Neural Information Processing Systems_ , volume 32, 2019.
  * Li et al. (2019) Gen Li, Xiaokang Liu, and Kun Chen.  Integrative multi-view regression: Bridging group-sparse and low-rank models.  _Biometrics_ , 75(2):593–602, 2019.
  * Loukas et al. (2021) Lefteris Loukas, Manos Fergadiotis, Ion Androutsopoulos, and Prodromos Malakasiotis.  EDGAR-CORPUS: Billions of tokens make the world go round.  In _Proceedings of the Third Workshop on Economics and Natural Language Processing_ , pages 13–18, Punta Cana, Dominican Republic, 2021.
  * Lounici et al. (2011) Karim Lounici, Massimiliano Pontil, Sara Van De Geer, and Alexandre B Tsybakov.  Oracle inequalities and optimal inference under group sparsity.  _The Annals of Statistics_ , 39(4):2164–2204, 2011.
  * Pilanci and Wainwright (2016) Mert Pilanci and Martin J Wainwright.  Iterative hessian sketch: Fast and accurate solution approximation for constrained least-squares.  _Journal of Machine Learning Research_ , 17(1):1842–1879, 2016.
  * Reinsel et al. (2022) Gregory C. Reinsel, Raja P. Velu, and Kun Chen.  _Multivariate Reduced-Rank Regression_.  Springer New York, NY, 2022.
  * Richtárik and Takáč (2016) Peter Richtárik and Martin Takáč.  Distributed coordinate descent method for learning with big data.  _Journal of Machine Learning Research_ , 17(75):1–25, 2016.
  * Ruhe (1970) Axel Ruhe.  Perturbation bounds for means of eigenvalues and invariant subspaces.  _BIT Numerical Mathematics_ , 10:343–354, 1970.
  * Salton and Buckley (1988) Gerard Salton and Christopher Buckley.  Term-weighting approaches in automatic text retrieval.  _Information Processing and Management_ , 24(5):513–523, 1988.
  * Temlyakov (2000) V. N. Temlyakov.  Weak greedy algorithms.  _Advances in Computational Mathematics_ , 12(2):213–227, 2000.
  * Temlyakov (2015) V. N. Temlyakov.  Greedy approximation in convex optimization.  _Constructive Approximation_ , 41(2):269–296, 2015.
  * Vershynin (2018) Roman Vershynin.  _High-Dimensional Probability: An Introduction with Applications in Data Science_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2018.
  * Wang et al. (2017) Jialei Wang, Jason D. Lee, Mehrdad Mahdavi, Mladen Kolar, and Nathan Srebro.  Sketching meets random projection in the dual: A provable recovery algorithm for big and high-dimensional data.  _Electronic Journal of Statistics_ , 11(2):4896–4944, 2017.
  * Wang et al. (2016) Xiangyu Wang, David Dunson, and Chenlei Leng.  Decorrelated feature space partitioning for distributed sparse regression.  In _Proceedings of the 30th International Conference on Neural Information Processing Systems_ , NIPS’16, pages 802–810, 2016.
  * Wedin (1972) Per-Åke Wedin.  Perturbation bounds in connection with singular value decomposition.  _BIT Numerical Mathematics_ , 12(1):99–111, 1972\.
  * Yang et al. (2016) Jiyan Yang, Michael W. Mahoney, Michael A. Saunders, and Yuekai Sun.  Feature-distributed sparse regression: A screen-and-clean approach.  In _Proceedings of the 30th International Conference on Neural Information Processing Systems_ , NIPS’16, pages 2711–2719, 2016.
  * Yeh et al. (2020) Hsiang-Yuan Yeh, Yu-Ching Yeh, and Da-Bai Shen.  Word vector models approach to text regression of financial risk prediction.  _Symmetry_ , 12(1), 2020.
  * Zheng et al. (2018) Wenjie Zheng, Aurélien Bellet, and Patrick Gallinari.  A distributed frank-wolfe framework for learning low-rank matrices with the trace norm.  _Machine Learning_ , 107(8):1457–1475, 2018.
  * Zhuo et al. (2020) Jiacheng Zhuo, Qi Lei, Alex Dimakis, and Constantine Caramanis.  Communication-efficient asynchronous stochastic frank-wolfe over nuclear-norm balls.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1464–1474, 2020.