# Class Symbolic Regression:
Gotta Fit ’Em All

Wassim Tenachi Université de Strasbourg, CNRS, Observatoire astronomique de
Strasbourg, UMR 7550, F-67000 Strasbourg, France  Rodrigo Ibata Université de
Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000
Strasbourg, France  Thibaut L. François Université de Strasbourg, CNRS,
Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France
Foivos I. Diakogiannis Data61, CSIRO, Kensington, WA 6155, Australia

(Received December 04 2023)

###### Abstract

We introduce “Class Symbolic Regression” a first framework for automatically
finding a single analytical functional form that accurately fits multiple
datasets - each governed by its own (possibly) unique set of fitting
parameters. This hierarchical framework leverages the common constraint that
all the members of a single class of physical phenomena follow a common
governing law. Our approach extends the capabilities of our earlier Physical
Symbolic Optimization ($\Phi$-SO) framework for Symbolic Regression, which
integrates dimensional analysis constraints and deep reinforcement learning
for symbolic analytical function discovery from data. We demonstrate the
efficacy of this novel approach by applying it to a panel of synthetic toy
case datasets and showcase its practical utility for astrophysics by
successfully extracting an analytic galaxy potential from a set of simulated
orbits approximating stellar streams.

Symbolic regression, Reinforcement learning, Recurrent neural network, Class
symbolic regression, Hierarchical symbolic regression

††journal: ApJL

##  1 Introduction

Figure 1: Class Symbolic Regression framework sketch: searching for a unique
functional form simultaneously fitting multiple datasets. The process starts
at the left hand side, a batch of trial class analytical expressions are
generated using our $\Phi$-SO framework (Tenachi et al., 2023). The free
parameters appearing in those expressions are then optimized in a dataset-
specific manner i.e. allowing each dataset to have its own unique associated
values for each parameter. The neural network used to generate the trial
expressions is then reinforced based on the fit quality of the trial symbolic
functions. This process is repeated until convergence.

Since the beginning of the scientific revolution, researchers have tried to
find repeatable regularities in experiments and observations. Mathematical
structures were used in this exploration, and many new ones including
functions and differential equations were developed to respond to this need to
model nature. Perhaps because of shared symmetries between nature and
mathematics, these abstract structures have often been found to work
exceedingly well in reproducing and predicting properties of the world, to the
point where some have even considered whether the universe is actually
mathematical at heart (Tegmark, 2008).

The Symbolic Regression (SR) that the present contribution is concerned with
has a long pedigree. Perhaps its most famous application was by Kepler to
planetary ephemerides, thereby finding the fitting law that bears his name
(Kepler, 1609). This empirical law gave the observational basis upon which
Newton was able to build the physical theories developed in his Principia
Mathematica (Newton, 1687).

Modern SR (Schmidt & Lipson, 2009, 2011; Kommenda et al., 2020; Kammerer et
al., 2020; Bartlett et al., 2023b; Brence et al., 2021; Jin et al., 2019; Luo
et al., 2022; Tohme et al., 2023; Udrescu & Tegmark, 2020; Udrescu et al.,
2020; Kamienny et al., 2022; Biggio et al., 2020, 2021; Vastl et al., 2022;
Kamienny et al., 2023; Martius & Lampert, 2017; Brunton et al., 2016; Zheng et
al., 2022; Sahoo et al., 2018; Petersen et al., 2021a; Landajuela et al.,
2022; Holt et al., 2023; Scholl et al., 2023) aims to use the immense
computational resources at our disposal to search through possible analytic
descriptions in terms of a set of functions and operators (e.g. $x$, $+$, $-$,
$\times$, $/$, $\sin$, $\cos$, $\exp$ $\log$, …) to best fit some numerical
dataset $(\mathbf{x},y)$ we wish to model. Concretely, one seeks some analytic
function $f:\mathbb{R}^{n}\longrightarrow\mathbb{R}$ that fits
$y=f(\mathbf{x})$ given those data. It is worth pointing out here that the
search space becomes exponentially larger the longer the analytic expression
is that we seek to find. Hence the key to SR is to develop efficient schemes
to search through the possibilities, and most importantly, to prune out poor
choices.

Our modern computational abilities have allowed us to examine nature in
unprecedented quantitative detail, with cameras, spectrographs and other
detectors amassing vast quantities of numerical data. It is likely that the
clues to next-generation physics and understanding lie therein, and so we are
tasked to devise methodologies capable of handling this wealth of information
and translating it into coherent, interpretable and intelligible physical
models. The promise of SR is that it may allow us in part to answer this need
to find accurate and intelligible empirical laws in giant datasets to best
capitalize on the community’s observational investments.

While SR has been extensively applied in scientific research, its focus has
largely been on single dataset analysis, overlooking the rich potential in
examining multiple datasets linked to a singular physical phenomenon. The
present article extends our $\Phi$-SO framework (presented in Tenachi et al.
2023) further by allowing the search for a function that can simultaneously
fit several datasets at once, each with (possibly) different fitting
parameters. This opens up the new possibility of implementing a functional
search on the properties of a _class_ of objects. This approach is relevant
across various natural sciences, but it particularly shines in astrophysics,
where multiple observations of a single phenomenon are often available,
providing a rich multi-dataset setup enabling us to devise ’universal’ laws
that apply to a range of celestial objects of interest.

The layout of the paper is as follows: In Section 2, we present the
methodology of our approach. In Section 3 we apply this method to synthetic
data from a panel of scientific cases. In Section 4, we illustrate the
practical application of our method in a more complex scenario by
demonstrating the recovery of a Milky Way-like analytic galactic potential,
which encapsulates the dark matter distribution, from simulated orbits that
could be inferred from stellar streams. Finally, Section 5, offers a
discussion and a conclusion.

##  2 Method

We build our “Class Symbolic Regression” framework on the Physical Symbolic
Optimization ($\Phi$-SO) framework for Symbolic Regression. This framework
combines deep reinforcement learning with in situ dimensional analysis
constraints to construct solutions that avoid physically nonsensical
combinations of units. This algorithm currently achieves state-of-the-art
performance on physics datasets, and significantly outperforms competitors on
the standard Feynman SR benchmark La Cava et al. (2021) in exact symbolic
recovery in the presence of even slight levels of noise (exceeding $0.1\%$).

Figure 1 gives an overview of our Class SR framework. Using $\Phi$-SO we
generate a batch of analytical expressions via a recurrent neural network
(RNN). Subsequently, we individually optimize the free parameters in these
expressions for each of the $N$ datasets $\\{(\mathbf{x},y)_{i}\\}_{i<N}$,
assigning unique values $\\{\theta_{i}\\}_{i<N}$ associated to each dataset.
This optimization is conducted using the L-BFGS nonlinear optimization routine
Zhu et al. (1997). Encoding our mathematical symbols with PyTorch Paszke et
al. (2019), enables us to use PyTorch’s implementation of the L-BFGS routine,
which benefits from PyTorch’s auto-differentiation capabilities to efficiently
optimize these free parameters.

We then use reinforcement learning to update the RNN’s parameters following a
risk-seeking gradient policy (Petersen et al., 2021a), as detailed in Tenachi
et al. (2023). This update is based on the average quality of the fit of the
trial functional form $f$ across datasets, evaluated using a normalized root
mean squared error (NRMSE):
$\text{NRMSE}=\frac{1}{\sigma_{y}}\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_{i}-f(\mathbf{x}_{i}))^{2}}$.
We repeat this process until the RNN converges to a unique high quality
expression and its associated parameter values simultaneously fit all
datasets.

It should be noted that it is straightforward for our framework to accommodate
any number of class-parameters - that is, parameters shared across the entire
class and having similar values across all datasets — in addition to the
dataset-specific parameters. This is achieved through a two step L-BFGS search
at each optimization iteration: a first step focuses on improving the class-
parameters, while second step is dedicated to fine-tuning the dataset-specific
parameters.

Furthermore, the sequential nature of expression generation in our $\Phi$-SO
framework enables the incorporation of various priors regarding the resulting
expressions as demonstrated in (Tenachi et al., 2023; Bartlett et al., 2023a;
Petersen et al., 2021b; Kim et al., 2021). This allows for customized
constraints on the generated expressions such as adherence to the rules of
dimensional analysis (which was one of the focal points of our previous study
Tenachi et al. 2023) but also simpler priors such as constraints on the number
of occurrences of given parameters, the length of the expression and more.

##  3 Multi-Dataset SR Challenges

Case |  | Target
---
Expression
# Expressions | # Datasets |  | # Samples
---
per dataset
Variables |  | Parameters
---
(one set per dataset)
Harmonic Oscillator | $A.\cos(f.t+\Phi)$ | 100k | 4 | 8 | $t\in[0,3\pi]$ |  | $f\in\\{1.06,1.57,0.24,0.76\\}$
---
$\Phi\in\\{0.97,0.92,0.88,1.13\\}$
$A\in\\{1.03,1.16,0.67,1.02\\}$
Radioactive Decay | $n_{0}.e^{-t/T}$ | 100k | 4 | 16 | $t\in[0.5,6]$ |  | $n_{0}\in\\{0.47,0.98,1.34,1.76\\}$
---
$T\in\\{1.35,0.87,1.12,0.97\\}$
Free Fall | $-\frac{1}{2}.g.t^{2}+v_{0}.t+z_{0}$ | 100k | 4 | 8 | $t\in[0,1]$ |  | $g=9.81$
---
$v_{0}\in\\{7.34,-1.17,5.74,2.12\\}$
$z_{0}\in\\{1.23,6.48,-2.13,1.42\\}$

Table 1: Class Symbolic Regression challenges. Each row details a distinct challenge, with the objective being the exact symbolic recovery of the designated target expression using multiple synthetic datasets. Each dataset is generated by sampling the specified number of points from the target expression within the given variable ranges, using unique parameter sets for each dataset. The limit on expression evaluations during the search process is indicated in the corresponding column. |  | Recovery rate
---|---|---
Expression | Noise level | Class-SR | SR
$A.\cos(f.t+\Phi)$ | 0% | $4/4$ | $4/4$
10% | $4/4$ | $3/4$
$n_{0}.e^{-t/T}$ | 0% | $4/4$ | $3/4$
5% | $4/4$ | $2/4$
$-\frac{1}{2}.g.t^{2}+v_{0}.t+z_{0}$ | 0% | $4/4$ | $3/4$
5% | $4/4$ | $0/4$
| Average | 100 % | 62.5 %
Table 2: Comparison of exact symbolic recovery rates between Class Symbolic
Regression and standard Symbolic Regression. This table demonstrates the
enhanced effectiveness of Class-SR in identifying common underlying functions
across multiple datasets with varying scale parameters, resulting in a higher
success rate compared to the traditional SR method exploiting only one dataset
at a time.

Despite existing research efforts to establish benchmarks for SR (La Cava et
al., 2021; Matsubara et al., 2022; Graham et al., 2013), a benchmark tailored
specifically for Class SR has yet to be developed, reflecting the innovative
nature of this approach. To address this, we introduce our own Class SR
challenges, designed to evaluate a system’s ability to analyze multiple
datasets. These datasets represent varied observations of a similar phenomenon
occurring at different scales but governed by a consistent functional form.
Table 1 outlines these challenges, each focusing on accurately recovering the
symbolic expression from synthetic datasets having varied scale parameter
values. To heighten the challenge, the last scenario incorporates a parameter
that is common across all datasets in addition to other dataset-specific
parameters.

We run our Class SR algorithm using the hyper-parameters detailed in Tenachi
et al. (2023) but with dimensional analysis disabled to make the challenge
easier for other algorithms. We allow the use of the following operations:
$\\{+,-,\times,/,1/{\square},\sqrt{{\square}},{\square}^{2},-{{\square}},\exp,\log,\cos,\sin\\}$
as well as three adjustable free constants $\theta=\\{a,b,c\\}$ allowed to
have dataset-specific values and a constant equal to one $\\{1\\}$ but no
class-parameters. The recovery rate is evaluated by examining each expression
in the Pareto front, which represents the best fitting expressions at each
level of complexity generated by our algorithm. Successful recovery is noted
if an expression on the Pareto front is symbolically equivalent to the target
expression. Exact symbolic recovery is assessed by comparing these expressions
with the target expression using the SymPy library for symbolic mathematics
Meurer et al. (2017), following a methodology similar to the one in the
SRBench standardized SR benchmarking protocol (La Cava et al., 2021).

Table 2 presents a comparison of exact symbolic recovery rates between our
Class Symbolic Regression framework and the traditional SR approach under both
noiseless and noisy conditions, with noise modeled as per SRBench:
$y_{\text{noise}}=y+\epsilon,\quad\epsilon\sim\mathcal{N}\left(0,\gamma\sqrt{\frac{1}{N}\sum_{i}y_{i}^{2}}\right)$
where $\gamma$ is the level of noise. The Class SR recovery rate is determined
by varying the seed in multiple runs, while the regular SR rate is calculated
by analyzing each dataset independently.

Our results demonstrate that utilizing multiple datasets generally enhances
the recovery rate, particularly in noisy scenarios where noise overfitting is
a concern. Notably, our Class SR algorithm effectively recovered the free fall
expression, including identifying that one of the dataset-specific parameters
actually remains constant across all datasets. While prior knowledge of the
presence of such class-parameters alongside dataset-specific ones is
beneficial as it constrains the search space, this showcases the robustness of
our algorithm, which can potentially identify class-parameters even without
predefined knowledge of their existence.

##  4 Recovering an analytic potential form stellar streams

We now turn to an astrophysical application of the algorithm: to try to find
the underlying potential of a gravitational system from a set of orbit
segments within it. This could be practically applicable for finding an
analytic potential model of a galaxy from a set of stellar streams. These
linear structures form from the tidal dissolution of globular clusters and
dwarf satellite galaxies. When their progenitors are of low mass, the escaping
stars have similar energy to the progenitor, and therefore follow a similar
orbit. Hence stellar streams approximate orbits in the host galaxy. As has
recently been shown by Ibata et al. (2023), for many real streams one can
calculate a “correction function” to convert an orbit model into a stream
track, and these functions are relatively insensitive to the adopted
potential. This procedure can be inverted to give the orbit from the stream.

For this test we imagine having access to full 6-dimensional phase-space
measurements of a sample of streams. For each structure $i$, the kinetic
energy per unit mass $E_{i,\rm kin}({\bf x})$ is simply:

$\frac{1}{2}{\bf v}^{2}=E_{\rm total}^{i}-\Phi({\bf x})\,.$ (1)

The total energy per unit mass $E_{\rm total}^{i}$, which is constant, but
different, for each stream, can be considered to be nuisance terms in our
search for the underlying potential $\Phi$. While it is possible to seek an
analytic form for $E_{i,\rm kin}({\bf x})$, by imposing solutions with a
constant dataset-specific additive offset by leveraging our $\Phi$-SO
framework’s capability to embed priors in the functional form or by allowing
only one dataset-specific parameter to model this offset, we opt for a more
challenging route. In this exercise, we intentionally avoid using this prior
knowledge. Instead, we deploy our algorithm by permitting only dataset-
specific parameters and exclude any class-parameters in order to showcase its
ability to deduce complex relationships without relying on predefined
assumptions.

For the present analysis we generated a sample of artificial orbit data that
approximates the sample of 29 thin and long streams studied by Ibata et al.
(2023). To this end we used the present day progenitor positions estimated by
Ibata et al. (2023), and integrated orbits within a universal (NFW) dark
matter halo model (Navarro et al., 1997) that very roughly approximates the
large-scale mass distribution in the Milky Way. The adopted potential (Łokas &
Mamon, 2001) is

$\Phi_{NFW}=-M_{200}\,g\ln[1+r/R_{h}](R_{h}/r)\,,$ (2)

where $M_{200}$ is the virial mass of the halo,
$g\equiv(\ln(1+c)-c/(1+c))^{-1}$ is a function of the halo concentration $c$
and $R_{h}$ is the scale radius. We chose $M_{200}=10^{12}{\rm\,M_{\odot}}$,
$c=10$ and $R_{h}=20.0{\rm\,kpc}$. The orbits consist of 100 phase space
points at locations between $\pm 1{\rm\,Gyr}$ from the current progenitor
location.

In this scenario, despite the presence of multiple alternative functional
forms that provide a near-perfect fit to individual streams, our Class SR
algorithm converges toward the correct functional form by exploring under 300k
expressions, correctly identifying an offset parameter specific to each stream
(corresponding to $E_{\rm total}^{i}$) and a functional form parameterized by
class-parameters common to all streams corresponding to $\Phi_{NFW}$,
showcasing its capabilities in practical astrophysical scenarios.

##  5 Discussion and conclusions

We presented a first framework for discovering symbolic analytical functions
that simultaneously fit multiple datasets by allowing for (possibly) unique
dataset-specific parameter values. This new framework which we dub “Class
Symbolic Regression” is built upon our earlier $\Phi$-SO framework which
already delivers state-of-the-art performances in symbolic recovery in the
presence of noise.

We demonstrated the efficacy of Class SR through various simple textbook
physics examples, finding better performance in exact symbolic recovery over
traditional SR, especially in noisy situations. Additionally, we applied our
method to a more complex astrophysical scenario, successfully rediscovering an
NFW galaxy potential model from orbits approximating stellar streams.

Regular SR, when applied to a single dataset, often risks overfitting to
specific characteristics of an observation, such as observational biases or
transient events, and noise. In contrast, our Class SR framework should
facilitate the finding of _universal_ analytical laws that apply to a range of
observations within a single class of physical phenomena. This enables our
framework to model the underlying physics rather than the specifics of
individual observations, with dataset-specific free parameters modeling the
unique aspects of each observation.

It should be noted that while Class SR might superficially resemble regular SR
applied to unbalanced datasets with dataset-specific parameters being akin to
additional input variables, this comparison is not entirely accurate. In Class
SR, these additional degrees of freedom represent unknown values that must be
determined, differentiating it as a distinct problem with its own unique
challenges.

A persistent issue in SR is model selection as the correct expression can
often be overlooked in favor of those that fit better or are less complex
(these concerns led to e.g., the development of single objective criterion
Bartlett et al. 2023b). Our framework, by searching for expressions that fit
multiple datasets, effectively utilizes information about the physical
phenomena’s class structure. This approach significantly mitigates model
selection challenges, helping avoid incorrect model choices influenced by
dataset-specific peculiarities. In addition, exploiting multiple datasets with
regular SR techniques would require fitting the individual datasets
independently, and then identifying the solutions in common between the
objects, which may not be possible if the measurements are uncertain.

In future work, we intend to improve on the machine learning aspects of our
method to more effectively leverage multiple datasets. As each dataset might
distinctly highlight certain symbolic terms or sub-expressions more
prominently than others, a promising strategy could be to periodically shift
the neural network’s training emphasis between datasets. This technique could
potentially refine the performance of Class SR by sequentially learning
different segments of the expression, rather than attempting to learn the
entire expression simultaneously, thereby facilitating the learning process.

## Code availability

The documented code for the $\Phi$-SO algorithm along with demonstration
notebooks are available on GitHub github.com/WassimTenachi/PhySO .

## Acknowledgments

RI acknowledges funding from the European Research Council (ERC) under the
European Unions Horizon 2020 research and innovation programme (grant
agreement No. 834148). The authors would like to acknowledge the High
Performance Computing Center of the University of Strasbourg for supporting
this work by providing scientific support and access to computing resources.
Part of the computing resources were funded by the Equipex Equip@Meso project
(Programme Investissements d’Avenir) and the CPER Alsacalcul/Big Data.

## References

  * Bartlett et al. (2023a) Bartlett, D., Desmond, H., & Ferreira, P. 2023a, in Proceedings of the Companion Conference on Genetic and Evolutionary Computation, GECCO ’23 Companion (New York, NY, USA: Association for Computing Machinery), 2402–2411, doi: 10.1145/3583133.3596327
  * Bartlett et al. (2023b) Bartlett, D. J., Desmond, H., & Ferreira, P. G. 2023b, IEEE Transactions on Evolutionary Computation, 1, doi: 10.1109/TEVC.2023.3280250
  * Biggio et al. (2020) Biggio, L., Bendinelli, T., Lucchi, A., & Parascandolo, G. 2020, in Learning Meets Combinatorial Algorithms at NeurIPS2020
  * Biggio et al. (2021) Biggio, L., Bendinelli, T., Neitz, A., Lucchi, A., & Parascandolo, G. 2021, in Proceedings of Machine Learning Research, Vol. 139, Proceedings of the 38th International Conference on Machine Learning, ed. M. Meila & T. Zhang (PMLR), 936–945.  https://proceedings.mlr.press/v139/biggio21a.html
  * Brence et al. (2021) Brence, J., Todorovski, L., & Džeroski, S. 2021, Knowledge-Based Systems, 224, 107077
  * Brunton et al. (2016) Brunton, S. L., Proctor, J. L., & Kutz, J. N. 2016, Proceedings of the National Academy of Sciences, 113, 3932, doi: 10.1073/pnas.1517384113
  * Graham et al. (2013) Graham, M. J., Djorgovski, S., Mahabal, A. A., Donalek, C., & Drake, A. J. 2013, Monthly Notices of the Royal Astronomical Society, 431, 2371
  * Holt et al. (2023) Holt, S., Qian, Z., & van der Schaar, M. 2023, in The Eleventh International Conference on Learning Representations.  https://openreview.net/forum?id=o7koEEMA1bR
  * Ibata et al. (2023) Ibata, R., Malhan, K., Tenachi, W., et al. 2023, arXiv preprint arXiv:2311.17202
  * Jin et al. (2019) Jin, Y., Fu, W., Kang, J., Guo, J., & Guo, J. 2019, arXiv preprint arXiv:1910.08892
  * Kamienny et al. (2023) Kamienny, P., Lample, G., Lamprier, S., & Virgolin, M. 2023, in Proceedings of Machine Learning Research, Vol. 202, International Conference on Machine Learning, ICML 2023, 23-29 July 2023, Honolulu, Hawaii, USA, ed. A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, & J. Scarlett (PMLR), 15655–15668.  https://proceedings.mlr.press/v202/kamienny23a.html
  * Kamienny et al. (2022) Kamienny, P.-A., d’Ascoli, S., Lample, G., & Charton, F. 2022, in Advances in Neural Information Processing Systems, ed. A. H. Oh, A. Agarwal, D. Belgrave, & K. Cho.  https://openreview.net/forum?id=GoOuIrDHG_Y
  * Kammerer et al. (2020) Kammerer, L., Kronberger, G., Burlacu, B., et al. 2020, in Genetic Programming Theory and Practice XVII (Springer), 79–99
  * Kepler (1609) Kepler, J. 1609, Astronomia nova.
  * Kim et al. (2021) Kim, J. T., Landajuela, M., & Petersen, B. K. 2021, in 1st Mathematical Reasoning in General Artificial Intelligence, International Conference on Learning Representations (ICLR)
  * Kommenda et al. (2020) Kommenda, M., Burlacu, B., Kronberger, G., & Affenzeller, M. 2020, Genetic Programming and Evolvable Machines, 21, 471
  * La Cava et al. (2021) La Cava, W., Orzechowski, P., Burlacu, B., et al. 2021, in Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks, ed. J. Vanschoren & S. Yeung, Vol. 1 (Curran).  https://datasets-benchmarks-proceedings.neurips.cc/paper_files/paper/2021/file/c0c7c76d30bd3dcaefc96f40275bdc0a-Paper-round1.pdf
  * Landajuela et al. (2022) Landajuela, M., Lee, C. S., Yang, J., et al. 2022, Advances in Neural Information Processing Systems, 35, 33985
  * Łokas & Mamon (2001) Łokas, E. L., & Mamon, G. A. 2001, MNRAS, 321, 155, doi: 10.1046/j.1365-8711.2001.04007.x
  * Luo et al. (2022) Luo, C., Chen, C., & Jiang, Z. 2022, International Journal of Computational Methods, 19, 2142002
  * Martius & Lampert (2017) Martius, G., & Lampert, C. H. 2017, Extrapolation and learning equations.  https://openreview.net/forum?id=BkgRp0FYe
  * Matsubara et al. (2022) Matsubara, Y., Chiba, N., Igarashi, R., & Ushiku, Y. 2022, in NeurIPS 2022 AI for Science: Progress and Promises.  https://openreview.net/forum?id=oKwyEqClqkb
  * Meurer et al. (2017) Meurer, A., Smith, C. P., Paprocki, M., et al. 2017, PeerJ Computer Science, 3, e103
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Newton (1687) Newton, I. 1687, Philosophiae Naturalis Principia Mathematica., doi: 10.3931/e-rara-440
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., et al. 2019, Advances in neural information processing systems, 32
  * Petersen et al. (2021a) Petersen, B. K., Larma, M. L., Mundhenk, T. N., et al. 2021a, in International Conference on Learning Representations.  https://openreview.net/forum?id=m5Qsh0kBQG
  * Petersen et al. (2021b) Petersen, B. K., Santiago, C., & Landajuela, M. 2021b, in 8th ICML Workshop on Automated Machine Learning (AutoML).  https://openreview.net/forum?id=yAis5yB9MQ
  * Sahoo et al. (2018) Sahoo, S., Lampert, C., & Martius, G. 2018, in International Conference on Machine Learning, PMLR, 4442–4450
  * Schmidt & Lipson (2009) Schmidt, M., & Lipson, H. 2009, science, 324, 81
  * Schmidt & Lipson (2011) —. 2011, Age-Fitness Pareto Optimization (New York, NY: Springer New York), 129–146, doi: 10.1007/978-1-4419-7747-2_8
  * Scholl et al. (2023) Scholl, P., Bieker, K., Hauger, H., & Kutyniok, G. 2023, arXiv preprint arXiv:2310.05537
  * Tegmark (2008) Tegmark, M. 2008, Foundations of Physics, 38, 101, doi: 10.1007/s10701-007-9186-9
  * Tenachi et al. (2023) Tenachi, W., Ibata, R., & Diakogiannis, F. I. 2023, arXiv e-prints, arXiv:2303.03192, doi: 10.48550/arXiv.2303.03192
  * Tohme et al. (2023) Tohme, T., Liu, D., & YOUCEF-TOUMI, K. 2023, Transactions on Machine Learning Research.  https://openreview.net/forum?id=lheUXtDNvP
  * Udrescu et al. (2020) Udrescu, S.-M., Tan, A., Feng, J., et al. 2020, Advances in Neural Information Processing Systems, 33, 4860
  * Udrescu & Tegmark (2020) Udrescu, S.-M., & Tegmark, M. 2020, Science Advances, 6, eaay2631
  * Vastl et al. (2022) Vastl, M., Kulhánek, J., Kubalík, J., Derner, E., & Babuška, R. 2022, arXiv preprint arXiv:2205.15764
  * Zheng et al. (2022) Zheng, W., Sharan, S., Fan, Z., et al. 2022, arXiv preprint arXiv:2212.14849
  * Zhu et al. (1997) Zhu, C., Byrd, R. H., Lu, P., & Nocedal, J. 1997, ACM Transactions on mathematical software (TOMS), 23, 550


# Impact of Physical Activity on Quality of Life During Pregnancy: A Causal ML
Approach

Kianoosh Kazemi1, Iina Ryhtä2, Iman Azimi3, Hannakaisa Niela-Vilén2, Anna
Axelin2, Amir M. Rahmani3, Pasi Liljeberg1
1Department of Computing, University of Turku, Finland,
2Department of Nursing Science, University of Turku, Turku, Finland
3Department of Computer Science, University of California, Irvine, USA
{kianoosh.k.kazemi, iikrry, hmniel, anmaax<EMAIL_ADDRESS>{azimii,
<EMAIL_ADDRESS>
###### Abstract

The concept of Quality of Life (QoL) refers to a holistic measurement of an
individual’s well-being, incorporating psychological and social aspects.
Pregnant women, especially those with obesity and stress, often experience
lower QoL. Physical activity (PA) has shown the potential to enhance the QoL.
However, pregnant women who are overweight and obese rarely meet the
recommended level of PA. Studies have investigated the relationship between PA
and QoL during pregnancy using correlation-based approaches. These methods aim
to discover spurious correlations between variables rather than causal
relationships. Besides, the existing methods mainly rely on physical activity
parameters and neglect the use of different factors such as maternal (medical)
history and context data, leading to biased estimates. Furthermore, the
estimations lack an understanding of mediators and counterfactual scenarios
that might affect them. In this paper, we investigate the causal relationship
between being physically active (treatment variable) and the QoL (outcome)
during pregnancy and postpartum. To estimate the causal effect, we develop a
Causal Machine Learning method, integrating causal discovery and causal
inference components. The data for our investigation is derived from a long-
term wearable-based health monitoring study focusing on overweight and obese
pregnant women. The causal graph is generated using the causal discovery
method and modified with the help of a domain expert team to accommodate the
mediators in the causal model. The machine learning (meta-learner) estimation
technique is used to estimate the causal effect. Our result shows that
performing adequate physical activity during pregnancy and postpartum improves
the QoL by units of 7.3 and 3.4 on average in physical health and
psychological domains, respectively. In the final step, four refutation
analysis techniques are employed to validate our estimation.

###### Index Terms:

Quality of Life, Physical Activity, Causal Machine Learning, Causal Discovery

##  I Introduction

Quality of Life (QoL) is a comprehensive concept that provides a holistic
evaluation of an individual’s well-being, emphasizing the biopsychosocial
aspects of health [1]. As defined by the World Health Organization (WHO), QoL
reflects an individual’s perspective on their life, taking into account
personal goals and expectations. Research indicates that pregnant women
generally experience lower QoL compared to population in general [2].

Factors such as obesity, stress, anxiety, and sleep disturbances can
negatively impact on QoL during pregnancy, while engaging in regular physical
activity (PA) can enhance it [2]. The significance of improving the QoL of
pregnant women lies in its potential to influence the well-being not only of
the mother but also of the entire family. It is recommended that pregnant
women aim to achieve 150 minutes of moderate PA per week. However, it has been
studied that [3] pregnant women with overweight and obesity met the PA
recommendations rarely. In this regard, enhancing the PA with overweight and
obese pregnant women is important since PA is known to have many short- and
long-term benefits for health, as well as for more holistic outcomes, e.g.,
QoL.

Several studies have been conducted to examine how PA impacts pregnant women’s
quality of life. For example, descriptive statistics measures, such as mean
value, standard deviation, confidence intervals, and standardized mean
differences, were used to explore the associations between PA and the QoL [4,
5, 6, 7]. Additionally, Pearson correlation coefficients have been utilized to
explore the association between physical activity and QoL [8]. In this study,
the Pregnancy Physical Activity Questionnaire tool was used to assess and
measure the physical activity levels of 141 pregnant women. In another study
[9], the relationship between PA and the QoL of pregnant women was assessed by
employing nonparametric tests such as one-way analysis of variance (ANOVA),
Chi-square test, and linear regression model.

There are several significant challenges involved in relying solely on
correlation analysis when analyzing a complex phenomenon with multiple causal
pathways. The causal relationships between variables must be disentangled
rather than being reliant on spurious correlations. In other words, the
conventional methods for estimating the effect of PA on QoL during pregnancy
rely mainly on correlation between variables and are inadequate to assess the
true underlying values. Moreover, most of the conventional techniques neglect
the use of different factors such as maternal (medical) history and context
data, including background information. Subsequently, these techniques
generate biased estimates. Many studies [10] have explored the limitations of
these methods, pointing out the need to go beyond just identifying data
associations and understanding the processes that generate this data and the
causal relationships between variables. Statistical approaches have their
limitations in revealing the direct and indirect causes linking pregnant
women’s QoL, PA, and other variables and in identifying and adjusting for
potential biases.

Furthermore, the existing studies relying on self-report questionnaires [6,
11] face limitations, including subjective reporting influenced by recall and
social desirability bias, as well as limited precision in responses.
Particularly in assessing pregnant women’s PA, subjective methods like
interviews and self-report questionnaires may inaccurately reflect PA. To
address these challenges, there is a growing trend toward the use of various
eHealth monitoring devices [12]. Wearables, measuring physical and health
conditions objectively [13, 14], provide continuous and accurate data,
overcoming the drawbacks of traditional subjective methods. This shift
enhances the quality of data, reduces dropout rates, and contributes to better
representation of health-related data on the levels and changes of activity
during pregnancy and the postpartum period. We believe that a Causal Machine
Learning (CML) method is required to exploit multi-modal data collection, such
as wearable-derived and subjective data (i.e., maternal (medical) history and
context data), while minimizing the bias of estimates.

In this paper, we examine the causal relationship between QoL and the PA of
pregnant women who are overweight and obese. To collect the required data, a
long-term health monitoring of 48 pregnant women was conducted. In this study,
we develop a CML approach consisting of two main components: causal discovery
and causal inference. In the first step, the causal discovery is exploited to
recover the underlying causal relationship between PA (the treatment variable)
and QoL (the outcome variable) utilizing the construction of causal graphs.
Then, we use the knowledge of the nursing expert team to modify the
constructed causal graph that clearly defines the set of variables pertinent
to both the PA and QoL. In the next step, the estimation algorithm (i.e., as a
part of the causal inference) is utilized to calculate the average treatment
effects for our dataset. By applying the CML method, the causal impact of PA
on the QoL during pregnancy and post-partum in this collected dataset is
examined. Our hypotheses are further validated through a variety of refutation
techniques, allowing us to draw well-informed conclusions.

##  II Method

In this section, we first present our data collection methodology, followed by
an explanation of the CML method. This study is part of a large intervention
study where lifestyle intervention was developed and tested with overweight
pregnant women [15]. In this paper, our focus is on the physical health and
psychological health domains of QoL.

###  II-A Data Collection

A long-term health monitoring study of overweight and obese pregnant women was
used for this study. The aim is to examine the relationship between PA and
QoL. Tracking vital signs, PA, and sleep was done using the OURA ring while
participants engaged in their normal daily routines.

Recruitment started in April 2021 and ended in May 2023. Pregnant women with
overweight and obesity (BMI$\mathrm{>}$25) were recruited before pregnancy
week 15. After explaining the study details, researchers provided OURA rings
to participants with usage instructions. Data collection occurred through an
online platform that we developed through which researchers could send
questionnaires and surveys to participants at 15 gestational weeks, 34
gestational weeks, and 12 weeks post-delivery, along with continuous OURA ring
monitoring from the recruitment until three months postpartum. Data from 48
pregnant women were used in our analysis. The participants’ backgrounds are
summarized in Table I.

Data collection was conducted using the OURA ring [16], a commercially
available wearable ring designed for monitoring a range of health parameters
and generating personalized daily health scores. The recorded daily data were
transmitted to smartphones and web applications through a Bluetooth
connection. In addition, women answered QoL questionnaires. QoL is measured by
the World Health Organization Quality of Life (WHOQOL–BREF) questionnaire [1].
WHOQOL-BREF is a 26-item assessment tool categorized into four domains:
physical health, psychological health, social relationships, and environment,
scoring between 0-100; higher scores indicate higher QoL. Moreover, a
questionnaire measuring depressive symptoms (EPDS) [17] was also completed at
each timepoint via the developed platform.

TABLE I: background information Characteristic |  | All participants
---
(n=48)
| Active group
---
(n=22)
| Low active group
---
(n=26)
Age, Mean (SD) | 29.3 (4.6) | 29.8 (5.2) | 29.2 (4.0)
| BMI, Mean (SD)
---
30.69 (4.74) | 29.47 (4.51) | 31.80 (4.76)
| Married or living
---
with a partner
– relationship status–
(yes) n (%),
Mean (SD)
47 (97.9) | 22 (100.0) | 24 (96.0)
| Employment status
---
n (%) Mean (SD)
$\bullet$ Working
$\bullet$ Unemployed
$\bullet$ Student
$\bullet$ Other
| 38 (79.2)
---
3 (6.3)
2 (4.2)
5 (10.4)
| 21 (95.5)
---
0 (0.0)
1 (4.5)
0 (0.0)
| 17 (68.0)
---
3 (12.0)
0 (0.0)
5 (20.0)
| EPDS at recruitment,
---
Mean (SD)
| 6.1 (4.1)
---
| 5.4 (3.5)
---
| 6.6 (4.6)
---
| Number of Children,
---
Mean (SD)
1.20 (1.6) | 1.3 (1.5) | 1.14 (1.6)

We separated participants into two distinct, non-overlapping subgroups: the
active group and the low active group. Utilizing the OURA ring, PA data (such
as activity type, intensity, duration, step counts, and the date/time of
activity) were collected daily. We used the medium-intensity activity
parameter (i.e., minutes per day that participants have performed activities
with an intensity level equivalent to walking) collected by OURA to calculate
the weekly duration of medium-intensity activity for each participant.
Subsequently, a threshold of 150 minutes per week, in alignment with the
recommended duration of moderate PA for pregnant women [3], was established.
This classification resulted in 22 women being placed in the active group,
with the remaining 26 pregnant women assigned to the low active group. Fig. 1
presents a visual representation of the weekly active time for the two groups.
Fig. 1 highlights that the active group exhibits a longer duration of medium
activity compared to the low active group throughout the phases of pregnancy
and postpartum. Furthermore, in the vicinity of the due date, the activity
duration converges, reaching a comparable level for both groups. In this
study, the following physical activity attributes are used for the analysis of
causal effect:

Figure 1: The activity duration for the Active vs. the Low active group.

  * •

Step counts: Total number of steps taken by an individual per day.

  * •

Average-MET: Shows the daily average MET when the user is awake. METs
(Metabolic Equivalents) are used to measure the intensity and energy
expenditure of different physical activities.

This research followed the ethical principles outlined in the Declaration of
Helsinki and the Finnish Medical Research Act (No 488/1999). It received
collective approval from the Ethics Committee of the Hospital District of
Southwest Finland and the maternity clinics (Protocol 113/1801/2020).
Participants received comprehensive oral and written information prior to
providing their written informed consent.

###  II-B Causal Discovery

In this study, we aim to uncover the underlying causal relationships between
QoL and PA during pregnancy by developing and applying causal discovery
algorithms. These algorithms will enable us to produce concrete causal graphs
that provide insight into how these relationships operate. Prior to employing
the causal discovery algorithms, it is essential to normalize our data,
ensuring its adherence to a standardized range. The objectives of the
normalization process are three-fold, including enhancing uniformity of data,
improving convergence during analysis and mitigating scale-related challenges.
Furthermore, standardization is carried out to calculate the covariance
matrix, which also serves as the basis for generating causal graphs.

In this study, we utilized three causal discovery algorithms: Peter-Clark (PC)
[18, 19, 20], Greedy Equivalence Search (GES) [20], and Greedy Interventional
Equivalence Search (GIES) [20]. These algorithms were selected based on the
characteristics of the dataset and computational resources available to us.
PC, a constrained-based causal discovery algorithm, initiates from a fully-
connected graph and identifies potential direct causal relations among
variables by employing conditional independence tests. Scalability and
capability are two main features of the PCs algorithm, making this method
suitable for working with large datasets. The GES and GIES algorithms build
causal graphs from an empty graph, then add or remove edges based on a score,
such as the Bayesian Information Criterion (BIC) [20]. GES algorithm can be
used to analyze causal relationships of both linear and non-linear types with
a wide range of variables. GIES [20], a derivative of GES, is designed to
reduce computational costs.

The causal graphs generated by the causal discovery component were analyzed
and the paths were adjusted to identify strongly contributing mediators. PA
affects QoL through a number of mediator variables. A team of nursing experts’
input was utilized iteratively in this process.

###  II-C Causal Inference

Causal inference attempts to estimate the causal effect of a treatment on a
particular outcome while adjusting the effect for certain variables that could
bias the estimate. In this study, a multi-mediator analysis [18] is utilized
to estimate the causal impact of the treatment variable (physical activity) on
the outcome (QoL score), taking into account the mediator variables identified
in the causal inference component. Causal effects can be estimated using a
variety of machine learning methods [18, 19, 20]. Our approach involves
constructing a meta-learner method (i.e., T-learner, a machine-learning
predictive model using different models for each treatment by splitting the
data) for the multi-mediator analysis, including both the Average Treatment
Effect (ATE) and mediator analysis. In ATE, average differences in outcomes
are calculated between treatment and control groups to estimate the causal
impact of a treatment variable [18, 19, 20]. An important technique in causal
inference is mediating effects [18], identifying how variables are influenced
by each other. An analysis of the mediate effect investigates how two
variables ($\mathrm{X}$ and $\mathrm{Y}$) interact through the mediating
influence of a third variable ($\mathrm{Z}$). The mediator analysis includes
three critical components to estimate the overall effect of $\mathrm{X}$ on
$\mathrm{Y}$: 1) the Natural Direct Effect (NDE) of $X$ on $\mathrm{Y}$, 2)
the Natural Indirect Effect (NIE) of $\mathrm{X}$ on $\mathrm{Y}$ via mediator
$\mathrm{Z}$, and 3) the Total Effect (TE) defined by the structural causal
model. TE consists of the same calculations as the ATE, including mediator
calculations as well.

NDE is defined in (1) as the expected change in $Y$ resulting from shifting
$\mathrm{X}$ from $\mathrm{x}$ to $x^{\prime}$ while keeping the mediating
variables at their previous values (i.e., $X=x$) before the shift from
$\mathrm{x}$to $x^{\prime}$.

$NDE_{x,x^{\prime}}(Y)=\sum\left[E\left(Y\mid x^{\prime},z\right)-E(Y\mid
x,z)\right]P(z\mid x)$ (1)

When $\mathrm{X}$ is maintained constant at $X=x$, and $\mathrm{Z}$ (for each
sample row) changes to the value it would have reached if $\mathrm{X}$ had
been set to $\mathrm{X}$ = $x^{\prime}$, NIE is the expected change in
$\mathrm{Y}$. Equation (2) can be used to calculate this counterfactual
concept.

$NIE_{x},{}_{x^{\prime}}(Y)=\sum_{z}E(Y\mid x,z)\left[P\left(z\mid
x^{\prime}\right)-P(z\mid x)\right]$ (2)

Equation (3) calculates mediators for the TE of $\mathrm{X}$ on $\mathrm{Y}$
in a Directed Acyclic Graph. Unlike linear systems, non-linear systems have no
additive effect. To be more precise, TE represents the difference between
direct and indirect effects (from $\mathrm{X}$ = $x^{\prime}$ to $\mathrm{X}$
= $\mathrm{x}$) of a reverse transition [18]. In summary, in the case of a
transition between $\mathrm{X}$ = $\mathrm{x}$ and $\mathrm{X}$ = $x^{\prime}$
(e.g., two groups of active and low active), the TE measures the predicted
change in $\mathrm{Y}$ due to the change in $\mathrm{X}$.

$TE_{x,x^{\prime}}(Y)=NDE_{x,x^{\prime}}(Y)-NIE_{x,x^{\prime}}(Y)$ (3)

###  II-D Refutation Analysis

Our estimations are tested by refuting them against unverified assumptions to
determine their sensitivity and robustness. As part of the refutation or
counterargument approach, noise is added to the common cause variable or the
treatment is replaced with a random variable in order to challenge the
estimates. Performing this analysis helps identify potential biases or
confounding factors that may affect the stability and robustness of the
estimate. The p-value less than 0.05 indicates the failure of the refutation
test, meaning the estimator and assumptions are invalid. Four refutation
methods are used in our experiments [18]. Placebo Treatment: In this test, an
independent random variable replaces the true treatment variable. It is
expected that the new estimation goes toward zero; otherwise, the assumption
is invalid. Data Subsets: This entails substituting the given dataset with a
randomly selected subset. The expected result is that the new estimation
remains unchanged; any deviation indicates an invalid assumption. Add Random
Common Cause: An independent random variable as a common cause is added to the
dataset. The new estimate is expected to remain unchanged; any changes signify
an invalid assumption. Unobserved common cause: This approach introduces a
supplemental dataset with a common cause between treatment and outcome. By
adjusting values to match known impact strengths, the robustness is assessed
via minimal estimator fluctuations. The expected result is that the new
estimation remains unchanged. In addition, this refutation does not return a
p-value.

##  III RESULTS

This section presents the findings of our CML approach using different
refutation techniques. To estimate causal effects using the proposed CML
method, we employed four stages: 1) causal graph construction (output of
causal discovery component), 2) the estimand identification using the causal
graph, 3) causal impact computation (output of causal inference component),
and 4) the estimate validation by manipulating assumptions (refutation
analysis).

The objective of this study is to determine the causal relationship between PA
(treatment variable) and QoL score (outcome) during pregnancy and post-partum,
ranging from 0 to 100. To achieve this, we applied the CML method to the
dataset (containing 48 pregnant women data) described in Section II-A. We
leveraged the causal discovery method to construct the causal graph. Then, by
exploiting the knowledge of our nursing expert team, the causal graph was
readjusted, as shown in Fig. 2.

In the Causal graph (Fig. 2), we included multiple mediators consisting of
steps, average-met from the OURA ring dataset, and the number of children,
work situation, and relationship status from the background information as
well as depression level from the EPDS questionnaire. We computed the ATE
using meta-learners (T-Learners) estimators. A light gradient boosting
regressor (max depth = 2 and min child samples = 60) was used in the meta-
learner to calculate the outcomes for each patient. As shown in Table II,
causal estimates are presented for each QoL domain score (i.e., physical
health and psychological). As shown in this table, the ATE estimates for all
three time-points were positive, meaning that the Active group had higher QoL
scores in physical health. For example, the estimated effect of being
physically active in week 15 was calculated as 10, showing that being active
causes a 10-unit increase in perception of the physical health score domain
when pregnant women are considered active. Similarly, in week 34 of pregnancy,
the ATE value was 4.41, meaning that the active group had a better physical
health score by 4 units. In addition, in post-partum, an ATE of 7.03 is
achieved, showing an increase of 3 units in physical health score compared to
week 34.

Figure 2: Causal graph.

Likewise, in the psychological domain, in all three timepoints, the ATE value
was positive, showing that active pregnant women have reported a higher
psychological score. For instance, during weeks 15 and 34, the ATE values were
2.81 and 2.67, respectively, meaning that active pregnant women have higher
psychological scores during early pregnancy and mid-pregnancy by roughly 2
units. Moreover, it is worth mentioning that the ATE values for both physical
health and psychological domains have a higher magnitude in the first
timepoint (i.e., week 15 of pregnancy) and the last timepoint (i.e., post-
partum week 12) compared to week 34. This variation in ATE values is due to
the fact that the perception of QoL is affected strongly by the progression of
pregnancy [7]. Since QoL is likely to be affected by multiple factors,
determining individual effects is difficult.

###  III-A Refutation Analysis

In this study, four refutation analysis methods were used to validate the
estimated results obtained by causal discovery. Table II summarizes the
refutation results. Our estimated causal effects were validated by all four
refutation tests. Based on the p-values (which are higher than 0.05), it is
evident that there are no significant causal effects in the refutation
scenarios. Moreover, we were able to confirm the robustness of our initial
estimation by considering the new estimate for the placebo treatment, which
was approximately zero, and the new estimates for the other methods remain
unchanged.

In summary, the proposed CML showed having more PA during pregnancy and post-
partum increases the QoL score in a causal way. Our causal graph incorporated
the direct and indirect effects of the PA parameters as well as the medical
and context information collected subjectively to minimize the bias of
estimates. Moreover, to have a better understanding and meaningful
interpretation of QoL scores, we classified QoL scores as (100–95) Near
perfect QoL, (95–85) Very good QoL, (85–70) Good QoL, (70—57.5) Moderately
QoL, and (57.5–40) Somewhat bad QoL [21]. As shown in the result, there is an
average 7.3-unit increase in each physical health score from low active to
active group, indicating changes in an individual’s QoL experience (for
example, from Near perfect to Very good). In future work, we extend the CML
method to personalized treatment effects and incorporate transportability into
the model for the transfer of causal effects to a new population.

##  IV CONCLUSION

This paper explored the causal effect of PA on QoL among pregnant women.
Maternal (medical) history and context data, as well as wearable-derived data,
were also considered crucial factors in our analysis. In this regard, a CML
method was developed, including both causal discovery and causal inference
components. Based on our findings, we can make an inference about how the PA
affects pregnant women’s QoL. Our result indicated that the PA led to an
improvement in the QoL during pregnancy and post-partum. As part of our causal
graph generation efforts, we incorporated domain expertise.

TABLE II: Validation result for proposed causal graph Method | Physical Health Score | Psychological Score
---|---|---
| Estimated
---
effect
| New
---
estimated
effect
p-value |  | Estimated
---
effect
| New
---
estimated
effect
p-value
Week 15
| Placebo
---
treatment
10.08 | 0.48 | 0.38 | 2.81 | 0.18 | 0.42
| Add random cause
---
10.08 | 9.88 | 0.38 | 2.81 | 2.98 | 0.38
Data subset | 10.08 | 9.95 | 0.32 | 2.81 | 2.40 | 0.29
| Unobserved
---
random cause
10.08 | 10.04 | — | 2.81 | 2.82 | —
| Week 34
| Placebo
---
treatment
4.41 | 0.19 | 0.36 | 2.67 | -0.02 | 0.48
| Add random cause
---
4.41 | 0.41 | 0.45 | 2.67 | 2.70 | 0.82
Data subset | 4.41 | 3.90 | 0.25 | 2.67 | 2.26 | 0.25
| Unobserved
---
random cause
4.41 | 0.47 | — | 2.67 | 2.65 | —
| Post-Partum Week 12
| Placebo
---
treatment
7.03 | 0.01 | 0.32 | 5.04 | 0.40 | 0.46
| Add random cause
---
7.03 | 7.02 | 0.49 | 5.04 | 4.98 | 0.43
Data subset | 7.03 | 7.36 | 0.47 | 5.04 | 5.53 | 0.19
| Unobserved
---
random cause
7.03 | 7.06 | — | 5.04 | 5.82 | —

## References

  * [1] World Health Organization. WHOQOL: Measuring Quality of Life, available ”https://www.who.int/tools/whoqol”, 2012.
  * [2] Lagadec, N., et al. “Factors influencing the quality of life of pregnant women:a systematic review”,BMC Pregnancy Childbirth, 18, 455, 2018.
  * [3] Bacchi, E., et al. “Physical activity patterns in normal-weight and overweight/obese pregnant women.” PLoS One, vol. 11, e0166254, 2016
  * [4] Seneviratne, S.N. et al., “Effects of antenatal exercise in overweight and obese pregnant women on maternal and perinatal outcomes: A randomised controlled trial”, BJOG Int. J. Obstet. Gynaecol. vol. 123, pp. 588–597, 2016.
  * [5] O’Connor, P.J., et al., “Effects of resistance training on fatigue-related domains of quality of life and mood during pregnancy: A randomized trial in pregnant women with back pain”, Psychosom. Med. e327, 2018.
  * [6] Rodríguez-Blanque, R., et al., “Water exercise and quality of life in pregnancy: A randomised clinical trial”, Int. J. Environ. Res. Public Health, vol. 17, pp. 1288, 2020.
  * [7] Nascimento, S. L., et al. “The effect of an antenatal physical exercise programme on maternal/perinatal outcomes and quality of life in overweight and obese pregnant women: a randomised clinical trial.” BJOG, no. 118, vol. 12, pp. 1455-1463, 2011.
  * [8] Mourady, Danielle, et al. “Associations between quality of life, physical activity, worry, depression and insomnia: A cross-sectional designed study in healthy pregnant women.” PloS one, no. 12, vol. 5, 2017.
  * [9] Bahadoran, Parvin, and Soheila Mohamadirizi. “Relationship between physical activity and quality of life in pregnant women.” Iran. J. Nurs. Midwifery Res., no. 20, pp. 282-6, 2015.
  * [10] Sanchez, P., et al. “Causal machine learning for healthcare and precision medicine.” Royal Society Open Science vol.9, pp. 220638, 2022.
  * [11] Krzepota, J., et al., “Relationships between physical activity and quality of life in pregnant women in the second and third trimester.” Int. J. Environ. Res. Public Health, vol. 12, pp. 2745, 2018.
  * [12] Niela-Vilén, H., et al. “Pregnant women’s daily patterns of well-being before and during the COVID-19 pandemic in Finland: Longitudinal monitoring through smartwatch technology.”PloS one, e0246494, 2021.
  * [13] Saarikko, J., et al. “Continuous 7-Month Internet of Things–Based Monitoring of Health Parameters of Pregnant and Postpartum Women: Prospective Observational Feasibility Study.” JMIR Form Res.,2020.
  * [14] Sarhaddi, Fatemeh, et al. ”Long-term IoT-based maternal monitoring: system design and evaluation.” Sensors vol.21, no. 7, pp. 2281, 2021.
  * [15] Saarikko, Johanna, et al. “Supporting lifestyle change in obese pregnant mothers through the wearable internet-of-things (SLIM)-intervention for overweight pregnant women: Study protocol for a quasi-experimental trial.” PloS one 18, no.1, pp. e0279696, 2023.
  * [16] OURA ring, Available online: https://ouraring.com/.
  * [17] Cox JL, et al., “Detection of postnatal depression. Development of the 10-item Edinburgh Postnatal Depression Scale”, Brit. J. Psychiat., 1987.
  * [18] Pearl, J. ”Causal inference in statistics: An overview.”, pp. 96-146, 2009.
  * [19] Pearl, J. ”The causal mediation formula—a guide to the assessment of pathways and mechanisms.”, Prevention Science, pp. 426-436, 2012.
  * [20] Glymour, C., et al. “Review of causal discovery methods based on graphical models”, Frontiers in genetics, vol.10, pp. 524, 2019.
  * [21] Hyland, M. E., et al., ”Development of a new type of global quality of life scale, and comparison of performance and preference for 12 global scales.” Quality of Life Research, vol. 5, pp. 469-480, 1996.


$a$ is good for mean convergence along $\ell$-term APs for all totally ergodic
$\ell$-step nilsystems. enumerate Furthermore, both properties are implied by
the following one: enumerate

  3. (iii)

$a$ is good for $\ell$-step irrational equidistribution.

###### Proof.

The implication $(i)\implies(ii)$ is trivial and the implication
$(iii)\implies(ii)$ follows from Proposition 7.3. So it remains to establish
the implication $(ii)\implies(i)$. The argument is similar to the one used to
establish the implication $(ii)\implies(i)$ in Theorem 8.2, so we omit the
details. ∎

###  8.3. Proof of Theorem 1.3

We restate and prove Theorem 1.3.

######  Theorem 8.4 (Criteria for convergence along $\ell$-term APs -
nilsystems).

Let $a\colon{\mathbb{N}}\to{\mathbb{Z}}$ be a sequence and
$\ell\in{\mathbb{N}}$. Then the following properties are equivalent:

  1. (i)

$a$ is good for mean convergence along $\ell$-term APs for all nilsystems.

  2. (ii)

$a$ is good for mean convergence along $\ell$-term APs for all $\ell$-step
nilsystems.

  3. (iii)

$a$ is good for pointwise convergence along $\ell$-term APs for all
$\ell$-step nilsystems.

  4. (iv)

$a$ is good for $\ell$-step equidistribution.

###### Proof.

The equivalence of the last three properties was established in Proposition
7.2 and the implication $(i)\implies(ii)$ is trivial.

So it remains to establish the implication $(ii)\implies(i)$. Since the
seminorm $\lvert\\!|\\!|\cdot|\\!|\\!\rvert_{s+1}$ is a norm on an $s$-step
nilsystem [18, Chapter 12, Theorem 17], all sequences are good for seminorm
control along $\ell$-term APs for every nilsystem. Taking into account this
property, the argument is identical to the one used to establish the
implication $(ii)\implies(i)$ of Theorem 1.1 (see previous subsection). ∎

###  8.4. Proof of Theorem 1.5.

We restate and prove Theorem 1.5.

######  Theorem 8.5 (Criteria for convergence of square averages).

Let $a,b\colon{\mathbb{N}}\to{\mathbb{Z}}$ be sequences such that $a,b,a+b$
are good for seminorm control for all ergodic systems. Then the following
properties are equivalent:

  1. (i)

For all systems $(X,{\mathcal{X}},\mu,T)$ and functions $f_{1},f_{2},f_{3}\in
L^{\infty}(\mu)$, we have

(58) $\lim_{N\to\infty}\frac{1}{N}\sum_{n=1}^{N}\,T^{a(n)}f_{1}\cdot
T^{b(n)}f_{2}\cdot
T^{a(n)+b(n)}f_{3}=\lim_{N\to\infty}\frac{1}{N^{2}}\sum_{r,s=1}^{N}\,T^{r}f_{1}\cdot
T^{s}f_{2}\cdot T^{r+s}f_{3}.$

  2. (ii)

Identity (58) holds for all ergodic $2$-step nilsystems.

Furthermore, property $(i)$ holds for all totally ergodic systems if and only
if property $(ii)$ holds for all totally ergodic $2$-step nilsystems.

###### Proof.

The implication $(i)\implies(ii)$ is obvious.

We establish the implication $(ii)\implies(i)$. Using an ergodic decomposition
argument we can assume that the system is ergodic. So lets fix an ergodic
system $(X,{\mathcal{X}},\mu,T)$ and suppose that the averages on the left
hand side of (58) are good for degree $3$ seminorm control for this system. We
claim that they are then good for degree $2$ seminorm control. To see this,
note first that our assumption implies that

$\lim_{N\to\infty}\Big{(}{\mathbb{E}}_{n\in[N]}\,T^{a(n)}f_{1}\cdot
T^{b(n)}f_{2}\cdot
T^{a(n)+b(n)}f_{3}-{\mathbb{E}}_{n\in[N]}\,T^{a(n)}\tilde{f}_{1}\cdot
T^{b(n)}\tilde{f}_{2}\cdot T^{a(n)+b(n)}\tilde{f}_{3}\Big{)}=0$

where the limit is taken in $L^{2}(\mu)$ and
$\tilde{f}_{j}:={\mathbb{E}}(f_{j}|{\mathcal{Z}}_{2})$ for $j=1,2,3$. So it
suffices to show that if $\lvert\\!|\\!|f_{j}|\\!|\\!\rvert_{2}=0$ for some
$j=1,2,3$, then

$\lim_{N\to\infty}{\mathbb{E}}_{n\in[N]}\,T^{a(n)}\tilde{f}_{1}\cdot
T^{b(n)}\tilde{f}_{2}\cdot T^{a(n)+b(n)}\tilde{f}_{3}=0$

where the limit is taken in $L^{2}(\mu)$. By Theorem 8.1 the factor
${\mathcal{Z}}_{2}$ is an inverse limit of ergodic nilsystems, so using an
approximation argument we can assume that the system $(X,{\mathcal{X}},\mu,T)$
is an ergodic $2$-step nilsystem. In this case, property $(ii)$ gives that the
last limit is equal to

$\lim_{N\to\infty}{\mathbb{E}}_{r,s\in[N]}\,T^{r}f_{1}\cdot T^{s}f_{2}\cdot
T^{r+s}f_{3}.$

These averages are good for degree $2$ seminorm control, which follows by a
simple application of the Gowers-Cauchy-Schwarz inequality (see for example
[18, Chapter 8, Theorem 13]), hence the claim follows.

Thus, we have established that degree $3$ seminorm control for the averages on
the left hand side of (58) implies degree $2$ seminorm control and identity
(58) holds. Recall also that by assumption the averages on the left hand side
of (58) are good for seminorm control for the system
$(X,{\mathcal{X}},\mu,T)$. Hence, Theorem 1.11 applies and gives that the
averages on the left hand side of (58) are good for degree $2$ seminorm
control for this system and identity (58) holds.

The equivalence in the case where we restrict to totally ergodic systems in
$(i)$ and $(ii)$ is proved in a similar fashion. We omit the details. ∎

###  8.5. Proof of Theorem 1.4

For the proof of Theorem 1.4 we need a few preparatory results.

######  Lemma 8.6.

Let $a\colon{\mathbb{N}}\to{\mathbb{Z}}$ be a sequence such that for every
$\ell$-step nilmanifold $Y=H/\Delta$ with $H$ connected and simply connected
and all equidistributed polynomial sequences $(g(n)\cdot e_{Y})$ in $Y$, the
sequence $(g(a(n))\cdot e_{Y})$ is equidistributed in $Y$. Then the sequence
$a\colon{\mathbb{N}}\to{\mathbb{Z}}$ is good for $\ell$-step irrational
equidistribution.

######  Remark.

The converse is also true and based on the fact that polynomial sequences can
be represented as linear sequences on a larger nilmanifold. We will not use
this, though, so we do not prove it.

###### Proof.

Let $X=G/\Gamma$ be a connected $\ell$-step nilmanifold. Let $G_{0}$ be the
connected component of $G$. Then, as explained in Section 7.1, we have
$G=G_{0}\Gamma$ and $X$ can be identified with the nilmanifold
$X_{0}:=G_{0}/\Gamma_{0}$ where $\Gamma_{0}:=G_{0}\cap\Gamma$ (see [18,
Chapter 11, Section 4.2]). Let $\pi\colon X\mapsto X_{0}$ denote an
isomorphism.

Let $b\in G$ be arbitrary. Since $G=G_{0}\Gamma$, as explained in Section 7.1
the sequence $b^{n}$ can be factorized as $b^{n}=g_{0}(n)\gamma(n)$ where
$g_{0}(n)$ and $\gamma(n)$ are polynomial sequences in $G_{0}$ and $\Gamma$
respectively. It follows that for every $F\in C(X)$, there exists $F_{0}\in
C(X_{0})$ (in fact $F_{0}:=F\circ\pi$) such that

(59) $\int F\,dm_{X}=\int F_{0}\,dm_{X_{0}}\quad\text{ and }\quad F(b^{n}\cdot
e_{X})=F_{0}(g_{0}(n)\cdot e_{X_{0}}),\quad n\in{\mathbb{N}},$

and conversely, for every $F_{0}\in C(X_{0})$ there exists $F\in C(X)$ (in
fact $F:=F_{0}\circ\pi^{-1}$) such that (59) holds. Suppose now that $b\in G$
is ergodic in $X$ in which case by unique ergodicity we have

$\lim_{N\to\infty}{\mathbb{E}}_{n\in[N]}\,F(b^{n}\cdot e_{X})=\int F\,dm_{X}$

for every $F\in C(X)$. Then (59) gives that

$\lim_{N\to\infty}{\mathbb{E}}_{n\in[N]}\,F_{0}(g_{0}(n)\cdot e_{X_{0}})=\int
F_{0}\,dm_{X_{0}}$

for every $F_{0}\in C(X_{0})$. Hence, the sequence $(g_{0}(n)\cdot e_{X_{0}})$
is equidistributed in $X_{0}$. Since $X_{0}=G_{0}/\Gamma_{0}$ where $G_{0}$ is
connected and simply connected, our assumption applies (for $Y:=X_{0}$,
$H:=G_{0}$, $\Delta:=\Gamma_{0}$) and gives that the sequence
$(g_{0}(a(n))\cdot e_{X_{0}})$ is equidistributed in $X_{0}$. We deduce from
this and (59) that the sequence $(b^{a(n)}\cdot e_{X})$ is equidistributed in
$X$. This proves that the sequence $a\colon{\mathbb{N}}\to{\mathbb{Z}}$ is
good for $\ell$-equidistribution and completes the proof. ∎

The only use of this result for our purposes is to prove the next lemma. The
advantage is that in the case of a connected and simply connected nilpotent
Lie group $G$ the map $g\mapsto g^{k}$ is surjective for every
$k\in{\mathbb{N}}$.

######  Lemma 8.7.

Suppose that the sequence $(da(n))$ is good for $\ell$-step irrational
equidistribution for some $d\in{\mathbb{N}}$. Then the sequence $(a(n))$ is
also good for $\ell$-step irrational equidistribution

###### Proof.

Let $X=G/\Gamma$ be an $\ell$-step nilmanifold where $G$ is connected and
simply connected and $g(n)$ be a polynomial sequence in $G$ such that
$(g(n)\cdot e_{X})$ is equidistributed in $X$. By Lemma 8.6, it suffices to
show that the sequence $(g(a(n))\cdot e_{X})$ is equidistributed in $X$. Note
that since $G$ is connected, it is divisible (see [18, Chapter 10, Corollary
9]), meaning that every element has a $k$-th root in $G$ for every
$k\in{\mathbb{N}}$. Hence, there exists a sequence $\tilde{g}$ in $G$ such
that $\tilde{g}(dn)=g(n)$ for all $n\in{\mathbb{N}}$. It is easy to see that
$\tilde{g}$ is a polynomial sequence: since $g$ can be written down as a
product of sequences of the form $a^{n^{i}}$, the sequence $\tilde{g}$ will be
a product of sequences of the form $\tilde{a}^{n^{i}}$ for
$\tilde{a}:=a^{1/d^{i}}$ a $d^{i}$-th root of $a$.

We claim that the sequence $(\tilde{g}(n)\cdot e_{X})$ is also equidistributed
in $X$. Indeed, by Leibman’s equidistribution criterion (see Section 7.1) it
suffices to verify this property when $G$ is a finite dimensional torus. By
Weyl’s equidistribution criterion, this amounts to showing that if a
$p\in{\mathbb{R}}[n]$ is a polynomial and the sequence $(p(dn))$ is
equidistributed in ${\mathbb{T}}$, then so is the sequence $(p(n))$. This is
clearly the case since the first polynomial has a non-constant irrational
coefficient if and only if the second one does. This proves the claim.

Since the polynomial sequence $(\tilde{g}(n)\cdot e_{X})$ is equidistributed
in $X$ and by assumption the sequence $(da(n))$ is good for irrational
equidistribution, Lemma 8.6 gives that the sequence $(\tilde{g}(da(n))\cdot
e_{X})=(g(a(n))\cdot e_{X})$ is equidistributed in $X$. This completes the
proof. ∎

The next lemma was proved in [13, Lemma 7.6].

######  Lemma 8.8.

Let $d\in{\mathbb{N}}$. Suppose that the sequence $(a(n))$ is good for
$\ell$-step irrational equidistribution and the set
$I_{d}:=\\{n\in{\mathbb{N}}\colon d|a(n)\\}$ has positive density. Then the
sequence $a_{d}\colon{\mathbb{N}}\to{\mathbb{Z}}$, defined by $(a(n))_{n\in
I_{d}}$, is also good for $\ell$-step irrational equidistribution.

Combining the previous two lemmas, we get the following.

######  Lemma 8.9.

Let $d\in{\mathbb{N}}$. Suppose that the sequence $(a(n))_{n\in{\mathbb{N}}}$
is good for $\ell$-step irrational equidistribution and the set
$I_{d}:=\\{n\in{\mathbb{N}}\colon d|a(n)\\}$ has positive density. Then the
sequence $a_{d}\colon{\mathbb{N}}\to{\mathbb{Z}}$, defined by $(a(n))_{n\in
I_{d}}$, is also good for $\ell$-step irrational equidistribution.

We also need a similar fact about seminorm control.

######  Lemma 8.10.

Let $d,\ell,s\in{\mathbb{N}}$, $k_{1},\ldots,k_{\ell}\in{\mathbb{Z}}$ be non-
zero and distinct, and $a\colon{\mathbb{N}}\to{\mathbb{Z}}$ be a sequence.
Suppose that the sequences $k_{1}a,\ldots,k_{\ell}a$ are good for degree $s$
seminorm control for all ergodic systems, and the set
$I_{d}:=\\{n\in{\mathbb{N}}\colon d|a(n)\\}$ has positive density. Let
$a_{d}\colon{\mathbb{N}}\to{\mathbb{Z}}$ be the sequence defined by
$(a(n)/d)_{n\in I_{d}}$. Then the sequences $k_{1}a_{d},\ldots,k_{\ell}a_{d}$
are good for degree $s+1$ seminorm control for all ergodic systems.

###### Proof.

The argument is similar to the one used in the proof of [10, Lemma 5.2].

An ergodic decomposition argument shows that the sequences
$k_{1}a,\ldots,k_{\ell}a$ are good for degree $s$ seminorm control for all
systems (not necessarily ergodic).

We consider the system $(X_{d},{\mathcal{X}}_{d},\mu_{d},T_{d})$ where

$X_{d}:=X\times\\{0,\ldots,d-1\\},\quad{\mathcal{X}}_{d}:={\mathcal{X}}\otimes{\mathcal{P}}(\\{0,\ldots,d-1\\}),\quad\mu_{d}:=\mu\times\nu_{d}$

with $\nu_{d}:=\frac{\delta_{0}+\cdots+\delta_{d-1}}{d}$, and for $x\in X$, we
let

$T_{d}(x,i):=\begin{cases}(x,i+1),\,&i=0,\ldots,d-2,\\\
(Tx,0),\;&i=d-1.\end{cases}$

The key property is that $T_{d}^{d}(x,i)=(Tx,i)$ for all $x\in X$ and
$i\in\\{0,\ldots,d-1\\}$, hence for $f\in L^{\infty}(\mu)$, we have

(60) $(T_{d}^{d}(f\otimes 1))(x,i)=(Tf)(x).$

Applying our good seminorm assumption for the product system

$(X_{d}\times
X_{d},{\mathcal{X}}_{d}\otimes{\mathcal{X}}_{d},\mu_{d}\times\mu_{d},T_{d}\times
T_{d}),$

we get that if $g_{1},\ldots,g_{\ell}\in L^{\infty}(\mu_{d})$ and
$\lvert\\!|\\!|g_{j}\otimes\overline{g_{j}}|\\!|\\!\rvert_{s,T_{d}\times
T_{d}}=0$ for some $j\in[\ell]$, then

(61) $\displaystyle\lim_{N\to\infty}{\mathbb{E}}_{n\in[N]}\,(T_{d}\times
T_{d})^{k_{1}a(n)}(g_{1}\otimes\overline{g_{1}})\cdots(T_{d}\times
T_{d})^{k_{\ell}a(n)}(g_{\ell}\otimes\overline{g_{\ell}})=0.$

Using the bound (26), Lemma 4.6(ii) for $w_{n}:={\bf 1}_{I_{d}}(n)$,
$n\in{\mathbb{N}}$, and the fact that $I_{d}$ has positive density, we deduce
that

(62) $\lim_{N\to\infty}{\mathbb{E}}_{n\in
I_{d,N}}\,T_{d}^{k_{1}a(n)}g_{1}\cdots T_{d}^{k_{\ell}a(n)}g_{\ell}=0$

in $L^{2}(\mu_{d})$ whenever $\lvert\\!|\\!|g_{j}|\\!|\\!\rvert_{s+1,T_{d}}=0$
for some $j\in[\ell]$, where for $N\in{\mathbb{N}}$ we let
$I_{d,N}:=I_{d}\cap[N]$.

We claim that for this $s\in{\mathbb{N}}$, the sequence $(a(n)/d)_{n\in
I_{d}}$ is good for degree $s+1$ seminorm control along $\ell$-term APs for
the system $(X,{\mathcal{X}},\mu,T)$. To see this, let
$f_{1},\ldots,f_{\ell}\in L^{\infty}(\mu)$ be such that
$\lvert\\!|\\!|f_{j}|\\!|\\!\rvert_{s+1,T}=0$ for some $j\in[\ell]$, say for
$j=j_{0}$. It suffices to show that

$\lim_{N\to\infty}{\mathbb{E}}_{n\in I_{d,N}}\,T^{k_{1}a(n)/d}f_{1}\cdots
T^{k_{\ell}a(n)/d}f_{\ell}=0$

in $L^{2}(\mu)$. Let $g_{j}\in L^{\infty}(\mu_{d})$ be defined by
$g_{j}:=f_{j}\otimes 1$ for each $j=1,\ldots,\ell$, and recall that
$\mu_{d}=\mu\times\nu_{d}$. By (60), it suffices to show that

$\lim_{N\to\infty}{\mathbb{E}}_{n\in I_{d,N}}\,T_{d}^{dk_{1}a(n)/d}g_{1}\cdots
T_{d}^{dk_{\ell}a(n)/d}g_{\ell}=0$

in $L^{2}(\mu_{d})$, or equivalently that

$\lim_{N\to\infty}{\mathbb{E}}_{n\in I_{d,N}}\,T_{d}^{k_{1}a(n)}g_{1}\cdots
T_{d}^{k_{\ell}a(n)}g_{\ell}=0.$

By (62), this holds as long as we establish that
$\lvert\\!|\\!|g_{j_{0}}|\\!|\\!\rvert_{s+1,T_{d}}=0$. To see this, note that
since $g_{j_{0}}=f_{j_{0}}\otimes 1$ and $\mu_{d}=\mu\times\nu_{d}$, our
assumption $\lvert\\!|\\!|f_{j_{0}}|\\!|\\!\rvert_{s+1,T}=0$ and (60) give
that
$\lvert\\!|\\!|g_{j_{0}}|\\!|\\!\rvert_{s+1,T_{d}^{d}}=\lvert\\!|\\!|f_{j_{0}}|\\!|\\!\rvert_{s+1,T}=0$.
The implication $\lvert\\!|\\!|g_{j_{0}}|\\!|\\!\rvert_{s+1,T_{d}}=0$ then
follows from (24). This completes the proof. ∎

We are now ready to prove Theorem 1.4, we restate it for convenience.

######  Theorem 8.11.

If the sequence $a\colon{\mathbb{N}}\to{\mathbb{Z}}$ is good for seminorm
control along $\ell$-term APs, is good for $\ell$-step irrational
equidistribution, and has good divisibility properties, then it is good for
recurrence along $\ell$-term APs.

###### Proof.

Let $(X,{\mathcal{X}},\mu,T)$ be a system and $f\in L^{\infty}(\mu)$ be such
that $f\geq 0$ and $\int f\,d\mu>0$. Let also
$k_{1},\ldots,k_{\ell}\in{\mathbb{Z}}$. Our goal is to show that for some
subset $I\subset{\mathbb{N}}$ with positive density, we have

(63) $\lim_{N\to\infty}{\mathbb{E}}_{n\in I_{N}}\int f\cdot
T^{k_{1}a(n)}f\cdots T^{k_{\ell}a(n)}f\,d\mu>0,$

where as before, for $N\in{\mathbb{N}}$ we let $I_{N}:=I\cap[N]$. We can
assume that $k_{1},\ldots,k_{\ell}$ are non-zero and distinct and the system
is ergodic.

Since the sequence $(a(n))$ is good for seminorm control along $\ell$-term
APs, and the set $I$ is assumed to have positive density, Lemma 4.6(ii) for
$w_{n}:={\bf 1}_{I}(n)$, $n\in{\mathbb{N}}$, and the bound (26) imply that it
is also good for seminorm control of the averages in (63). Hence, using
Theorem 8.1, we can assume that the system is an inverse limit of ergodic
nilsystems, and using an approximation argument (as in [6, Lemma 2.1] or [16,
Lemma 3.2]), we can assume that it is an ergodic nilsystem. As explained in
Section 7.1, there exists $d\in{\mathbb{N}}$ such that $X$ can be partitioned
as $X=\bigcup_{i=0}^{d-1}T^{i}X_{0}$, and for $i=0,\ldots,d-1$, the
transformation $T^{d}$ leaves the set $T^{i}X_{0}$ invariant and the
corresponding system is totally ergodic.

Since $\int f\,dm_{X}>0$ and $X=\bigcup_{i=0}^{d-1}T^{i}X_{0}$, there exists
$i_{0}\in\\{0,\ldots,d-1\\}$ such that $\int f\,dm_{T^{i_{0}}X_{0}}>0$. We let

$I_{d}:=\\{n\colon d|a(n)\\}.$

We claim that for $I:=I_{d}$ the positivity property (63) holds, or
equivalently, that

(64) $\lim_{N\to\infty}{\mathbb{E}}_{n\in I_{d,N}}\int f\cdot
T^{k_{1}a(n)}f\cdots T^{k_{\ell}a(n)}f\,dm_{X}>0$

where for $N\in{\mathbb{N}}$ we let $I_{d,N}:=I_{d}\cap[N]$.

Let $S:=T^{d}$. For $i\in[d]$, the transformation $S$ leaves the set
$T^{i}X_{0}$ invariant, and the system $(T^{i}X_{0},m_{T^{i}X_{0}},S)$ is
totally ergodic. Let $(a_{d}(n))$ be the sequence defined by $(a(n)/d)_{n\in
I_{d}}$. Since the set $I_{d}$ has positive density, it suffices to show that

(65) $\lim_{N\to\infty}{\mathbb{E}}_{n\in[N]}\int f\cdot
S^{k_{1}a_{d}(n)}f\cdots S^{k_{\ell}a_{d}(n)}f\,dm_{T^{i_{0}}X_{0}}>0.$

Our assumptions and Lemma 8.9 and 8.10 imply that the sequence $(a_{d}(n))$ is
good for $\ell$-step irrational equidistribution and good for seminorm control
along $\ell$-term APs for all ergodic systems. So we are now in position to
use the implication $(iii)\implies(i)$ in Theorem 8.3. It gives that the
sequence $(a_{d}(n))$ is good for $\ell$-convergence for all totally ergodic
systems. We apply this for the totally ergodic nilsystem
$(T^{i_{0}}X_{0},m_{T^{i_{0}}X_{0}},S)$, and we get that the average in (65)
equals

$\lim_{N\to\infty}{\mathbb{E}}_{n\in[N]}\int f\cdot S^{k_{1}n}f\cdots
S^{k_{\ell}n}f\,dm_{T^{i_{0}}X_{0}}.$

Since $f\geq 0$ and $\int f\,dm_{T^{i_{0}}X_{0}}>0$, the last limit is
positive by Furstenberg’s multiple recurrence result [15]. This completes the
proof. ∎

## References

  * [1] M. Akcoglu, L. Sucheston, On operator convergence in Hilbert space and in Lebesgue space. Period. Math. Hungar. 2 (1972), 235–244.
  * [2] D. Berend, V. Bergelson. Mixing sequences in Hilbert spaces. Proc. Amer. Math. Soc. 98 (1986), 239–246
  * [3] V. Bergelson, B. Host, B. Kra, with an appendix by I. Ruzsa. Multiple recurrence and nilsequences. Inventiones Math. 160 (2005), no. 2, 261–303.
  * [4] V. Bergelson, A. Leibman. Polynomial extensions of van der Waerden’s and Szemerédi’s theorems. J. Amer. Math. Soc. 9 (1996), 725–753.
  * [5] V. Bergelson, A. Leibman. Cubic averages and large intersections. Recent trends in ergodic theory and dynamical systems. Contemporary Mathematics 631 (2015), 5–19.
  * [6] V. Bergelson, A. Leibman, E. Lesigne. Intersective polynomials and the polynomial Szemerédi theorem. Adv. Math. 219 (2008), no. 1, 369–388.
  * [7] N. Frantzikinakis. Multiple ergodic averages for three polynomials and applications. Trans. Amer. Math. Soc. 360 (2008), no. 10, 5435–5475.
  * [8] N. Frantzikinakis. Multiple correlation sequences and nilsequences. Inventiones Math. 202 (2015), 875–892.
  * [9] N. Frantzikinakis. Some open problems on multiple ergodic averages. Bulletin of the Hellenic Mathematical Society 60, (2016), 41–90.
  * [10] N. Frantzikinakis. Joint ergodicity of sequences. Preprint. arXiv:2102.09967
  * [11] N. Frantzikinakis, M. Johnson, E. Lesigne, M. Wierdl. Powers of sequences and convergence. Ergodic Theory Dynam. Systems 30 (2010), no. 5, 1431–1456.
  * [12] N. Frantzikinakis, B. Kuca. Joint ergodicity for commuting transformations and applications to polynomial sequences. Preprint. arXiv:2207.12288
  * [13] N. Frantzikinakis, E. Lesigne, M. Wierdl. Sets of $k$-recurrence but not $(k+1)$-recurrence. Ann. Inst. Fourier (Grenoble) 56 (2006), no. 4, 839–849.
  * [14] N. Frantzikinakis, E. Lesigne, M. Wierdl. Powers of sequences and recurrence. Proc. Lond. Math. Soc. (3) 98 (2009), no. 2, 504–530.
  * [15] H. Furstenberg. Ergodic behavior of diagonal measures and a theorem of Szemerédi on arithmetic progressions. J. Anal. Math. 71 (1977), 204–256.
  * [16] H. Furstenberg, Y. Katznelson. An ergodic Szemerédi theorem for commuting transformations. J. Analyse Math. 34 (1979), 275–291.
  * [17] B. Host, B. Kra. Non-conventional ergodic averages and nilmanifolds. Ann. of Math. 161 (2005), 397–488.
  * [18] B. Host, B. Kra. Nilpotent Structures in Ergodic Theory. Mathematical Surveys and Monographs, vol. 236. American Mathematical Society, Providence, RI, 2018.
  * [19] B. Host, B. Kra, A. Maass. Complexity of nilsystems and systems lacking nilfactors. J. Anal. Math. 124 (2014), 261–295.
  * [20] B. Host, B. Kra, A. Maass. Variations on topological recurrence. Monatsh. Math. 179 (2015), Issue 1, 57–89.
  * [21] T. Kamae, M. Mendès-France. Van der Corput’s difference theorem. Isr. J. Math. 31 (1978), 335–342.
  * [22] A. Leibman. Pointwise convergence of ergodic averages for polynomial sequences of rotations of a nilmanifold. Ergodic Theory Dynam. Systems 25 (2005), no. 1, 201–213.
  * [23] A. Leibman. Host-Kra and Ziegler factors and convergence of multiple averages. Handbook of Dynamical Systems, Vol. 1B, B. Hasselblatt and A. Katok, eds., Elsevier, (2006), 841–853.
  * [24] J. Leng. A quantitative bound for Szemerédi’s theorem for a complexity one polynomial progression over ${\mathbb{Z}}/N{\mathbb{Z}}$. Preprint (2022). arXiv:2205.05540
  * [25] J. Moreira, F. Richter. A spectral refinement of the Bergelson-Host-Kra decomposition and new multiple ergodic theorem. Ergodic Theory Dynam. Systems 39 (2019), no. 4, 1042-1070.
  * [26] S. Peluse. On the polynomial Szemerédi theorem in finite fields. Duke Math. J. 168 (2019), 749–774.
  * [27] S. Peluse. Bounds for sets with no polynomial progressions. Forum of Mathematics, Pi, 8 (2020), e16.
  * [28] S. Peluse, S. Prendiville. Quantitative bounds in the non-linear Roth theorem. Preprint 2019, arXiv:1903.02592
  * [29] A. Stepin. Flows on solvable manifolds. Upekhi Mat. Nauk 24 (1969), 241–242 (Russian).
  * [30] E. Szemerédi. On sets of integers containing no $k$ elements in arithmetic progression. Acta Arith. 27 (1975), 299–345.
  * [31] T. Ziegler. A nonconventional ergodic theorem for a nilsystem. Ergodic Theory Dynam. Systems 25 (2005), no. 4, 1357–1370.
  * [32] T. Ziegler. Universal characteristic factors and Furstenberg averages. J. Amer. Math. Soc. 20 (2007), 53–97.