\section{Introduction}

Dark matter (DM) is a key ingredient for models
that try to explain
cosmological structure formation without modifying gravity. Although the gravitational effects of DM are well documented \cite{Betoule:2014frx,Ade:2013zuv},
direct detection of particles for this dominant matter component continues to elude proofs: in accelerators or in nuclear recoil experiments
\cite{
Felcini:2018osp,
Klasen:2015uma},
indirect WIMP annihilation searches \cite{Conrad:2014tla}, in DM stars \cite{Dai:2009ik,Kouvaris:2015rea} or in some other indirect quests as illustrated in Refs.~\cite{
Bertolami:2012yp,
Abdalla:2007rd,
Delliou:2014awa}.

In this context, different testing avenues of possible DM effects are welcome, such as in pulsars, i.e. rotating neutron stars (NSs), which provide the advantage of extreme densities and can accrete DM, thus straining  the saturated neutron gas. The amount of DM acquired by a NS follows the Tolman-Oppenheimer-Volkoff (TOV) equation
, as in e.g. \cite{Tolos2015}.
Moreover, the effect of DM on NSs can directly lead to bounds for the masses of the different DM candidates \cite{Goldman:1989nd,Kouvaris:2011fi}.

Self-annihilating DM can also produce characteristic effects on NS
\cite{Kouvaris2008,
Bertone2008,
Kouvaris2010,
McCullough:2010ai,
deLavallaz2010,
PerezGarcia:2011hh,
PerezGarcia:2010ap}.
In particular, WIMPs annihilation in DM cores should produce temperature and luminosity changes, through heat, of old stars \cite{Kouvaris2008,Bertone2008,Kouvaris2010,deLavallaz2010}.
However, those changes  are
difficult to
detect \cite{Kouvaris2008,Sandin:2008db}.

In this short report, based on the results obtained in our previous work \cite{Deliyergiyev:2019vti} and on our recent paper \cite{DelPopolo:2019nng}, we propose a testable galactic probe for DM existence in the form of evolution of the pulsar mass towards the galactic centre (GC). According to the discussion above, NSs in increasingly DM rich environments should accrete more DM and thus display a characteristic mass decrease, the closer they are to the galactic centre.
We use NSs because first of all the very large baryon density inside NSs makes the interaction between baryons and DM following DM capture most likely; second, the NSs strong gravitational force makes DM particles escape very unlikely, after they interact and loose energy. This mass evolution is easier to test than other probes such as NS temperature time evolution with DM accretion, as discussed previously.

\section{Mass Change of Neutron Stars}
\label{sec:Implementation}

\begin{figure}[t]
  \begin{center}
    \includegraphics[width=0.395\textwidth]{densProfiles_Fig1a}
    \includegraphics[width=0.395\textwidth]{accreteMasses_Fig1b_3rdParam}
 \caption{
Figures taken from Ref.~\cite{DelPopolo:2019nng}.
(Left panel) Einasto profile for $\alpha=0.06,0.09,0.11,0.15,0.17$ (black dot dot dashed, blue dotted, red solid, yellow long dashed, and green dot dashed lines, respectively). The cyan dashed line is the Di Cintio profile \cite{DiCintio:2014xia}  (labelled DC14 in the text). (Right panel) The accreted mass according to
Kouvaris formula (Eq.~(\ref{eq:Kouv2013})). In this plot we do not consider the accreted mass corresponding to the Di Cinto profile, nor to
the Einasto profile with $\alpha=0.06$, in order to be conservative in our calculations.
}
\label{fig:einasto}
  \end{center}
\end{figure}

As shown by Ref.~\cite{Kouvaris2013}, in the framework of a spherically symmetric accretion scenario for a typical NS of mass $1.4 M_\odot$ and $R=10$ km,
the total accreted mass is given by Ref.~\cite{DelPopolo:2019nng}
\begin{equation}
M_{\rm acc}= 1.3 \times 10^{43} \left( \frac{\rho_{\rm dm}}{\rm 0.3 \, GeV/cm^3} \right) \left(\frac{\rm t}
{\rm Gyr} \right) f \,\,\, \rm GeV,
\label{eq:Kouv2013}
\end{equation}
which is an underestimation of a factor $\simeq 10$, since the accretion during the NS progenitor phase, of the same order as in the NS phase \cite{Kouvaris2010} (factor of 2), and the accretion coming from DM self-interaction \cite{Guver2014} are not taken into account.
Moreover, $M_{acc}$ of Eq.~(\ref{eq:Kouv2013}) should be corrected by factor 2.05 for 2 $M_{\odot}$

We obtain DM accretion $\simeq 10^{-11} {\rm M}_{\odot}$ for a typical NS in the solar neighborhood,
which is in agreement with the capture rates of Refs~\cite{Kouvaris2013, Guver2014, Zheng2016, Zhong2012} and the results from  \cite{Kouvaris2011}, but below the estimates from the DM accumulated using TOV.
A better agreement between the accreted DM mass and the accumulated DM  mass coming from TOV
is obtained for NSs located in Superdense DM clumps, Ultra Compact mini-haloes \cite{Berezinsky2013}, and close to the GC.

The another approach starts, as in usual structure formation, from linear perturbation, followed by non-linear collapse of DM, and then continued by baryons collapse in the DM potential wells previously formed \cite{Chang2018}.

Once the accreted DM mass, $M_{acc}$, as a function of the distance to the GC is determined from the right panel of Fig.~\ref{fig:einasto}, we can determine the corresponding mass change of the NS towards the galactic center.
In order to do so, we use Fig.10 of our previous work \cite{Deliyergiyev:2019vti}, where the maximum mass of NSs was obtained as a function of the DM mass inside the NS, $M_{DM}$, for the DM weakly interacting case, $y=0.1$\footnote{ The interaction strength is expressed in terms of the ratio of the DM fermion mass $m_f$, and scale of interaction $m_I$, $y=m_f/m_I$. This can be converted to usual units: one can estimate the cross section of DM self-interaction, taking the mass of the DM particle $m_f$ in units of GeV, as
\begin{align}
\sigma =& \frac{1}{4\pi}\frac{m_f^2}{m_I^4}
= \frac{y^4}{m_f^2} 3.2\times 10^{-27} \textrm{cm}^2  &
\rightarrow \sigma/m_f =& \frac{y^4}{m_f^3} 1.8\times 10^{-3} \frac{\textrm{cm}^2}{g}
\label{eq:yCrossSection}
\end{align}
}.
Moreover, we obtained similar plots for larger values of the interaction parameter $y$, $y=1,10,100$. Since $M_{DM}$ must be equal to the accreted mass, $M_{acc}$, we find a relation between the total mass of the NS, $M_{NS}$ ($M_T$ in the notation of \cite{Deliyergiyev:2019vti}) and the distance from the GC. The result is plotted in Fig.~\ref{fig:accreteMass}.
Its leftmost panel displays the change in mass of a NS moving towards the Milky Way (MW) halo center (colour and line coding as in Fig.~\ref{fig:einasto}), for a particle mass, $m_{DM}=500$ GeV, and $y=0.1$.
Our reference model ($\alpha=0.11$, red solid line) shows a NS mass change from 2 to 1 $M_{\odot}$ at 0.4 pc, while in the case $\alpha = 0.9$ (blue dotted line) the same change is observed at 1.35 pc. The other two cases show a slower mass change.
The yellow long dashed line ($\alpha = 0.15$), and the green dot dashed line ($\alpha = 0.17$) show that the mass reduces from 2 to 1.2 $M_{\odot}$, and from 2 to 1.7 $M_{\odot}$, at $10^{3}$ pc, respectively.
The next two panels show how the mass changes with decreasing the DM particle mass ($m_{DM}=200$ GeV, centre left, 100 GeV, centre right). Finally the rightmost plot shows the effect of the interaction strength for our reference profile ($\alpha=0.11$), and from right to left, in the cases $y=0.1;1;10;10;100$. Strong interaction ($y=100$) produces very small mass changes, while the weaker the interaction, the larger the change is.
\begin{figure}[t]
  \begin{center}
    \includegraphics[width=0.95\textwidth]{accreteMass_Fig2c_newlines}
    \caption{
Figures taken from Ref.~\cite{DelPopolo:2019nng}.
The changes of the NS mass from DM accumulation as a function of the NS distance to GC, for a particle mass $m_{DM}=500$ (leftmost), 200 (center left) and 100 (center right) GeV, with $y=0.1$ and  $\alpha=0.09; 0.11; 0.15; 0.17$ (blue dotted, red solid, yellow long dashed, and green dot dashed). The rightmost panel shows the role of the interaction strength for particle mass $m_{DM}=500$ GeV, $\alpha=0.11$, and $y=0.1, 1, 10, 100$ (solid, long dashed, dash dotted, and long dashed dotted lines).
    }
   \label{fig:accreteMass}
  \end{center}
\end{figure}

All the above results are obtained with conservative assumptions, i.e.
(a) not taking into account
(i) NS progenitor accretion, expected of order of the NS phase \cite{Kouvaris2010},
(ii) DM self-interaction accretion \cite{Guver2014},
(b) using a density profile from DM-only simulations, shallower than in recent hydrodynamic simulations \cite{DiCintio:2014xia}.

As the right panel of Fig.~\ref{fig:einasto} shows, we excluded from the analysis the $\alpha = 0.06$ Einasto profile, and even the DC14 profile \cite{DiCintio:2014xia}, a realistic profile calculated with hydro-dynamical simulations, to remain conservative. Although we made this choice, the orbital dynamics of PSR B1257+12 \cite{Iorio:2010hb}, together with the accretion predictions \cite{Bertone:2004pz} allows much larger DM accumulation in NS than in our present work: up to 10\% of a NS mass, in agreement with DC14 \cite{DiCintio:2014xia}. Moreover, complex astrophysical phenomena are occurring on sub-parsec scale near the GC, such as DM particles gravitational scattering by stars and capture in the supermassive BH, together with highly enhanced central density from the supermassive BH formation \cite{Gondolo1999}. Thus, more accurate density profiles can be introduced \cite{Bertone:2004pz}.

\section{Conclusion}

We have shown that the NS mass should reflect the changes of the DM environment in the MW(Fig.~\ref{fig:accreteMass}). This is done taking into account the DM accretion of NSs [see \cite{Deliyergiyev:2019vti}, for details], as it changes because of the increase of DM content \cite{Kouvaris2013} when we move towards the GC \cite{Bernal:2011pz}. This allows us to propose that the evolution of the pulsar
masses towards the GC of the MW can be a probe of the existence of DM. In fact, the decrease of the NS mass for NSs located closer and closer to the GC would put
constraints on the characteristics of the Galaxy halo dark matter profile, on the dark matter particle mass, and on the self-interaction strength. Such changes are expected to be observed in the near future in telescopes such as ngVLA \cite{murphy2018science}, SKA \cite{Konar:2016lgc}, Athena \cite{barcons2012athena, EESA2014athena}, NICER \cite{10.1117/12.926396} or eXTP \cite{Watts2018}.

\section{Acknowledgements}
M.D. work was supported by the Polish National Science Centre (NCN) grant 2016/23/B/ ST2 / 00692. MLeD acknowledges the financial support by Lanzhou University starting fund and the Fundamental Research Funds for the Central Universities (Grant No.lzujbky-2019-25). L.T. acknowledges support from the FPA2016-81114-P Grant from Ministerio de Ciencia, Innovacion y Universidades, Heisenberg Programme of the Deutsche Forschungsgemeinschaft under the Project Nr. 383452331 and PHAROS COST Action CA16214.


\section{Introduction}
Let $A \in \mathbb{R}^{m \times n}$ with $m \geq n$ and full rank. Suppose $x \in \mathbb{R}^n$ and
$$ Ax =b,$$
how would one go about finding $x$?
A popular method for solving these types of problems is the Kaczmarz method \cite{kac} (also known as the Algebraic Reconstruction Technique (ART) in computer tomography \cite{gordon, her, natterer} or as the Projection onto Convex Sets method \cite{cenker, sezan}). Denoting the rows of $A$ by $a_1, \dots, a_m$, the linear system can be written as a collection of inner products
$$ \forall~1 \leq i \leq m: \left\langle a_i, x\right\rangle = b_i.$$
The Kaczmarz method is an iterative scheme: given $x_k$, let us pick an equation, say the $i-$th equation, and modify $x_k$ the smallest possible amount necessary to make it correct: set $x_{k+1} = x_k + \delta a_i,$ where $\delta$ is such that $\left\langle a_i,x_{k+1} \right\rangle = b_i$. Formally,
$$ x_{k+1} = x_k - \frac{b_i - \left\langle a_i, x_k\right\rangle}{\|a_i\|^2}a_i.$$
The Kaczmarz method is then given by cycling through the indices. The method has a geometric interpretation (see Fig. 1): $x_{k+1}$ is obtained by projecting $x_k$ onto the hyperplane given by $\left\langle a_i, x \right\rangle = b_i$. This interpretation is already given in the original 1937 paper of Kaczmarz who ends his paper by saying that `The convergence of this method is completely obvious from a geometric perspective' \cite{kac}. However, given the intricate underlying geometry, convergence rates are difficult to obtain. In a seminal paper, Strohmer \& Vershynin \cite{strohmer} showed that by choosing the equations in a random order (where the $i-$th equation is chosen with a likelihood proportional to $\|a_i\|^2$), the method converges exponentially: more precisely, they show that
$$ \mathbb{E} \left\| x_k - x \right\|_2^2 \leq \left(1 - \frac{1}{\|A\|_F^2 \| A^{-1}\|_2^2}\right)^k \|x_0 - x\|_2^2,$$
where $\|A^{-1}\|_2$ is the operator norm of the inverse and $\|A\|_F$ is the Frobenius norm.
This is a remarkable result in many ways. The analysis of Kacmarz and Kacmarz-type algorithms has become quite popular: we refer to
\cite{bai, eldar, gower, lee, leventhal, liu, ma, need, need2, need3, need4, tan, zouz} and references therein.

\begin{center}
\begin{figure}[h!]
\begin{tikzpicture}[scale=2]
\draw[thick] (-1, -1) -- (1, 1);
\draw[thick] (0, -1) -- (0, 1);
\draw[thick] (1, -1) -- (-1, 1);
\filldraw (0,0) circle (0.03cm);
\node at (0.2, 0) {$x$};
\filldraw (0.7, 0.7) circle (0.03cm);
\node at (0.9, 0.7) {$x_k$};
\filldraw (0, 0.7) circle (0.03cm);
\node at (0.25, 0.78) {$x_{k+1}$};
\draw [ <-] (0.1, 0.7) -- (0.7, 0.7);
\filldraw (-0.7/2, 0.7/2) circle (0.03cm);
\filldraw (-0.35, 0.35) circle (0.03cm);
\node at (-0.6, 0.35) {$x_{k+2}$};
\draw [->] (0, 0.7) -- (-0.3, 0.4);
\end{tikzpicture}
\caption{A geometric interpretation: iterative projections onto the hyperplanes given by $\left\langle a_i, x\right\rangle = b_i$.}
\end{figure}
\end{center}

\section{ Results}
\subsection{Summary}
A priori, one might be inclined to believe that the very nature of the algorithm produces an essentially random sequence in the sense that
the normalized vectors $(x_k-x)/\|x_k - x\|$ behave mostly randomly
or, alternatively put, that the approximations converge to $x_k \rightarrow x$ `from all directions'. It is the goal of this short note to show
that this is decidedly \textit{not} the case: the singular vectors $v_1, \dots, v_n$ of the matrix $A$ describe directions of distinguished dynamics. We describing three main observations:
\begin{enumerate}
\item \S 2.1: while $\|x_k-x\|$ decreases exponentially, this is happening at different rates in different subspaces. We show (Theorem 1) that for any (right) singular vector $v_{\ell}$ of $A$ associated to the singular value $\sigma_{\ell}$
$$\mathbb{E}  \left\langle x_{k}-x, v_{\ell}\right\rangle = \left(1 - \frac{\sigma_{\ell}^2 }{\|A\|_F^2}\right)^k\left\langle x_0-x, v_{\ell}\right\rangle.$$
Contributions from the smallest singular vectors ($\sigma_{\ell}$ small) decays the slowest: hence we expect them to be dominant for $k$
large. We recall that for the smallest singular vector $\sigma_n$, we have the identity
$$ \sigma_n = \|A^{-1}\|_2 \qquad \mbox{leading to the rate} \qquad 1 - \frac{1}{\|A\|_F^2 \|A^{-1}\|_2^2}$$
which recovers the convergence rate of Strohmer \& Vershynin \cite{strohmer}. Thus in the generic setting one
cannot hope to improve the rate.
\item \S 2.2: this phenomenon also impacts convergence rates:
if $x_k$ is not mainly the linear combination of singular vectors corresponding to small singular values, Randomized Kaczmarz
converges \textit{faster}: we have
$$\mathbb{E} ~ \|x_{k+1} - x\|_2^2 \leq \left(1 - \frac{1}{\|A\|_F^2} \left\| A  \frac{x_k-x}{\|x_k - x\|} \right\|^2 \right) \|x_{k} - x\|_2^2.$$
We observe that this factor may actually be quite large since
$$ \|A^{-1}\|^2\leq  \left\| A  \frac{x_k-x}{\|x_k - x\|} \right\|^2 \leq \|A\|^2.$$
The lower bound is only close to being attained when $x_k-x$ is close to the space spanned by singular vectors given by small singular values of $A$.
\item \S 2.3: once $x_k$ is mainly the combination of singular vectors corresponding to small singular values, the Randomized Kaczmarz Method is also less likely to explore adjacent subspaces. We show that the normalized vector $(x_k-x)/\|x_k-x\|$ then moves only very slowly on $\mathbb{S}^{n-1}$ since
$$ \mathbb{E} ~\left\langle \frac{x_k-x}{\|x_k-x\|}, \frac{x_{k+1}-x}{\|x_{k+1}-x\|}  \right\rangle^2 = 1 - \frac{1}{\|A\|_F^2}\left\| A  \frac{x_k-x}{\|x_k - x\|} \right\|^2.$$
\end{enumerate}

\subsection{Transfer to Small Singular Vectors} We start by establishing the main phenomenon which can be phrased as a simple algebraic identity.
\begin{theorem} Let $v_{\ell}$ be a (right) singular vector of $A$ associated to the singular value $\sigma_{\ell}$. Then
$$\mathbb{E} \left\langle x_{k} - x, v_{\ell} \right\rangle = \left(1 - \frac{\sigma_{\ell}^2 }{\|A\|_F^2}\right)^k\left\langle x_0 - x, v_{\ell}\right\rangle.$$
\end{theorem}
We recall that
$$ \sum_{i=1}^{n}{\sigma_{i}^2} = \|A\|_{F}^2,$$
so the arising factor multiplicative factor is always in between 0 and 1. This shows that the Randomized Kaczmarz method has exponential decay at different
exponential rates that depend on the singular values. In particular, the slowest rate of decay is given by the smallest singular value $\sigma_n$. Since
$$ \sigma_n = \|A^{-1}\|_2 \qquad \mbox{the rate is} \qquad 1 - \frac{1}{\|A\|_F^2 \|A^{-1}\|_2^2}$$
which recovers the convergence rate of Strohmer \& Vershynin \cite{strohmer}. In particular, if $x_0$ has contributions from singular vectors corresponding to large singular values, the sequence initially decays faster than expected since components corresponding to larger singular values decay quicker.
 A previous result in that spirit was obtained by Jiao, Jin \& Lu \cite{jiao} who obtained an `averaged' version of Theorem 1: by decomposing $\mathbb{R}^n = L + H$ into subspaces spanned by low singular vectors $L$ (corresponding to small singular values) and high singular vectors $H$ (corresponding to large singular values), they obtain the estimate
$$ \mathbb{E} \| \pi_{H} x_{k+1} \|^2 \leq (1-c_1) \| \pi_H x_{k}\|^2 + c_2 \| \pi_{L} x_k\|^2$$
where $\pi$ is the projection and $c_1, c_2 > 0$ are given in terms of singular values.

\subsection{Speed of Convergence}
The phenomenon described in Theorem 1 obviously has an effect on the speed of the algorithm: whenever $x_k - x$ is not close to singular vectors of $A$ that correspond to small singular vector, the Randomized Kaczmarz method actually converges faster. Theorem 2 is already implicit in Strohmer \& Vershynin \cite{strohmer} where a formulation of it appears as an intermediate step.
\begin{theorem} We have
$$\mathbb{E} ~ \|x_{k+1} - x\|_2^2 \leq \left(1 - \frac{1}{\|A\|_F^2} \left\| A  \frac{x_k-x}{\|x_k - x\|} \right\|^2 \right) \|x_{k} - x\|_2^2.$$
\end{theorem}
Observe that for any vector $\|v\|=1$, we have
$$ \|Av\|^2 \geq \sigma_n^2 \|v\|^2 = \|A^{-1}\|_2^2 \cdot  \|v\|_2,$$
where $\sigma_n$ is the smallest singular value. But clearly this last inequality is only sharp or close to sharp if $v$ is close to linear combination of singular vectors corresponding to small singular values: otherwise, the quantity is indeed larger.

\subsection{Stability of Small Singular Vectors.}
We conclude with a result of a geometric flavor. A natural object is the sequence of vectors
$$ \frac{x_k - x}{\|x_k - x\|} \in \mathbb{S}^{n-1}.$$
It essentially describes from which direction the sequence $(x_k)_{k=0}^{\infty}$ approaches the solution $x$. The result in the prior sections have shown that we expect this vector will eventually be concentrated on subspaces corresponding to singular vectors associated to small singular values. A different questions would be how much we expect the normalized vector to `move around'.
\begin{theorem}
If $x_k \neq x$ and $\mathbb{P}(x_{k+1} = x)=0$, then
$$ \mathbb{E} ~\left\langle \frac{x_k-x}{\|x_k-x\|}, \frac{x_{k+1}-x}{\|x_{k+1}-x\|}  \right\rangle^2 = 1 - \frac{1}{\|A\|_F^2}\left\| A  \frac{x_k-x}{\|x_k - x\|} \right\|^2.$$
\end{theorem}
We recall that for two unit vectors $\|u\| = 1 = \|v\|$, the inner quantity $\left\langle u, v\right\rangle^2$ can be interpreted as a notion of distance between the vectors: it is only close to 1 if $u \sim \pm v$. Since we only care about subspaces, the ambiguity in the sign does not play a role. We observe that the direction $(x_k-x)/\|x_k-x\|$ undergoes large fluctuations when $\|A (x_k-x)\|$ is large and is quite stable when $\|A(x_k-x)\|$ is small. Once $x_k$ is mainly comprised of small singular vectors, its direction does not change very much anymore: randomized Kaczmarz spends more time close to these subspaces than it does in other regions of space: the smaller the smallest singular value, the more pronounced phenomenon. Could this property be used for convergence acceleration?

\subsection{An Application: Finding Small Rayleigh Quotients}
We discuss an amusing application: finding vectors $\|x\|=1$ such that $\|Ax\|$ is small. In light of our results, such a sequence of vectors can be obtained from the Randomized Kaczmarz method applied to the equation
$$ Ax = 0.$$
We illustrate this with a simple example: we take $A \in \mathbb{R}^{1000 \times 1000}$ by picking each entry to be iid $\sim \mathcal{N}(0,1)$, we add $100\cdot \mbox{Id}_{1000}$ and then set the last row $a_{1000}$ to be a tiny perturbation of $a_{999}$ by adding 0.01 to each entry. This generates a matrix having singular
values between $\sim 0.5$ and $\sim 1.5$ and one tiny singular value at scale $\sim 10^{-4}$. We start with $x_0 = (1,\dots, 1)$ and plot $\|A x_k\|/\|x_k\|$ (see Fig. 2).
We observe that within $5000 \sim 5n$ iterations, the ratio is actually smaller than any singular value but the last one -- this implies that $x_k$ is actually already mostly given by $v_{1000}$. After 10.000 iterations, we observe $|\left\langle x_{10^4}/\|x_{10^4}\|, v_{1000}\right\rangle| \sim 0.9$.
\begin{center}
\begin{figure}[h!]
\begin{tikzpicture}[scale=2]
\node at (0,0.3) {\includegraphics[width=0.6\textwidth]{pic4.pdf}};
\end{tikzpicture}
\caption{A sample evolution of $\|Ax_k\|/\|x_k\|$.}
\end{figure}
\end{center}

\section{Proofs}
\subsection{Reduction to $Ax = 0$} We start with the ansatz
$$ x_k = x + y_k.$$
Plugging in, we obtain
\begin{align*}
x + y_{k+1} &= x_{k+1} = x_k + \frac{ b_i - \left\langle a_i, x_k\right\rangle}{\|a_i\|^2} a_i \\
&=x + y_k + \frac{ b_i - \left\langle a_i, x + y_k\right\rangle}{\|a_i\|^2} a_i \\
&= x + y_k  - \frac{  \left\langle a_i,  y_k\right\rangle}{\|a_i\|^2} a_i + \left( \frac{ b_i - \left\langle a_i, x\right\rangle}{\|a_i\|^2} a_i\right).
\end{align*}
We note that the final term in parenthesis vanishes because $x$ is the exact solution and thus $b_i - \left\langle a_i, x\right\rangle=0$. This results in
$$ y_{k+1} = y_k - \frac{  \left\langle a_i,  y_k\right\rangle}{\|a_i\|^2} a_i$$
which is exactly the Randomized Kaczmarz method for $Ay = 0$. Linearity of all the underlying objects (or, alternatively, the geometric interpretation in terms of hyperplane intersections) explains why this is perhaps not too surprising. In any case, this allows us to restrict ourselves to the case $Ax = 0$.
\subsection{Proof of Theorem 1}
\begin{proof}
As discussed above, we can assume without loss of generality that $b=0$. It remains to understand the behavior of the Randomized Kaczmarz method when applied to $Ax=0$. We have
\begin{align*}
\mathbb{E} \left\langle x_{k+1} , v_{\ell} \right\rangle &=
    \sum_{i=1}^{m}  \frac{\|a_i\|^2}{\|A\|_F^2} \left\langle    x_k - \frac{\left\langle a_i, x_k \right\rangle }{\|a_i\|^2} a_i  , v_{\ell} \right\rangle \\
   &=   \frac{1}{\|A\|_F^2}\sum_{i=1}^{m} \left( \|a_i\|^2 \left\langle x_k, v_{\ell}\right\rangle - \left\langle a_i, x_k\right\rangle \left\langle a_i, v_{\ell} \right\rangle \right)\\
   &= \left\langle x_k, v_{\ell}\right\rangle -   \frac{1}{\|A\|_F^2}\sum_{i=1}^{m} \ \left\langle a_i, x_k\right\rangle \left\langle a_i, v_{\ell} \right\rangle.
   \end{align*}
This sum has a special structure: $\left\langle a_i, x_k\right\rangle$ is simply the $i-$th entry of $A x_k$ and, likewise, $\left\langle a_i, v_{\ell}\right\rangle$ is simply the $i-$th entry of $A v_{\ell}$. This simplifies the sum. We then expand $x_k$ into right-singular vectors $v_i$, make use of $A v_i = \sigma_i u_i$ and the orthogonality of the left singular vectors $u_i$ to obtain
   \begin{align*}
 \left\langle x_k, v_{\ell}\right\rangle -   \frac{1}{\|A\|_F^2}\sum_{i=1}^{m} \ \left\langle a_i, x_k\right\rangle \left\langle a_i, v_{\ell} \right\rangle   &= \left\langle x_k, v_{\ell}\right\rangle - \frac{1}{\|A\|_F^2} \left\langle A x_k, A v_{\ell}\right\rangle \\
      &= \left\langle x_k, v_{\ell}\right\rangle- \frac{1}{\|A\|_F^2} \left\langle A \left(\sum_{i=1}^{n}{\left\langle x_k, v_i \right\rangle v_i}\right), A v_{\ell}\right\rangle \\
      &= \left\langle x_k, v_{\ell}\right\rangle- \frac{1}{\|A\|_F^2} \left\langle  \sum_{i=1}^{n}{\sigma_i \left\langle x_k, v_i \right\rangle u_i}, \sigma_{\ell} u_{\ell}\right\rangle\\&= \left\langle x_k, v_{\ell}\right\rangle\left(1 - \frac{\sigma_{\ell}^2 }{\|A\|_F^2}\right)
      \end{align*}
 which is the desired result.
\end{proof}
\subsection{Proof of Theorem 2}
\begin{proof} This proof follows the argument of Strohmer \& Vershynin \cite{strohmer}. If $Z$ is the (vector-valued) random variable assuming the value
$Z = a_j/\|a_j\|^2$ with probability $\|a_j\|^2/\|A\|_F^2$,
then Strohmer \& Vershynin have established the inequality
$$ \|x_{k+1} - x\|_2^2 \leq \left(1 - \left| \left\langle \frac{x_k-x}{\|x_k - x\|}, Z \right\rangle\right|^2 \right) \|x_{k} - x\|_2^2.$$
At this point we simply compute
\begin{align*}
\mathbb{E}  ~\left| \left\langle \frac{x_k-x}{\|x_k - x\|}, Z \right\rangle\right|^2 &= \sum_{i=1}^{m} \frac{\|a_j\|_2^2}{\|A\|_F^2} \left\langle \frac{x_k-x}{\|x_k - x\|},  \frac{a_j}{\|a_j\|_2} \right\rangle^2 \\
&=  \frac{1}{\|A\|_F^2} \sum_{i=1}^{m}   \left\langle \frac{x_k-x}{\|x_k - x\|},  a_j \right\rangle^2  \\
&= \frac{1}{\|A\|_F^2} \left\| A  \frac{x_k-x}{\|x_k - x\|} \right\|^2.
\end{align*}
\end{proof}

\subsection{Proof of Theorem 3}
\begin{proof} As above, we can assume without loss of generality that $b=0$. Moreover,
since the entire expression is invariant under scaling, we can assume w.l.o.g. that $\|x_k\|=1$. If $\mathbb{P}(x_{k+1}=0)=0$, we have
\begin{align*}
 \mathbb{E} \left\langle x_k, \frac{x_{k+1}}{\|x_{k+1}\|} \right\rangle^2 &= \sum_{i=1}^{m}  \frac{\|a_i\|^2}{\|A\|_F^2} \left\langle x_k,
   \frac{x_k - \frac{\left\langle a_i, x_k \right\rangle }{\|a_i\|^2} a_i}{\| x_k - \frac{ \left\langle a_i, x_k \right\rangle }{\|a_i\|^2} a_i \|}  \right\rangle^2 \\
   &= \sum_{i=1}^{m}  \frac{\|a_i\|^2}{\|A\|_F^2} \frac{ \left\langle x_k,
   x_k - \frac{\left\langle a_i, x_k \right\rangle }{\|a_i\|^2} a_i  \right\rangle^2}{\| x_k - \frac{ \left\langle a_i, x_k \right\rangle}{\|a_i\|^2} a_i \|^2}.
   \end{align*}
   We observe that
   $$x_k - \frac{ \left\langle a_i, x_k \right\rangle}{\|a_i\|^2} a_i  $$
   is the projection of $x_k$ onto the subspace $\left\langle x, a_i \right\rangle =0$ and thus
   $$ \left\langle x_k, x_k - \frac{ \left\langle a_i, x_k \right\rangle}{\|a_i\|^2} a_i  \right\rangle = \left\|x_k - \frac{ \left\langle a_i, x_k \right\rangle}{\|a_i\|^2} a_i  \right\|^2.$$
   Thus
   \begin{align*}
 \sum_{i=1}^{m}  \frac{\|a_i\|^2}{\|A\|_F^2} \frac{ \left\langle x_k,
   x_k - \frac{\left\langle a_i, x_k \right\rangle }{\|a_i\|^2} a_i  \right\rangle^2}{\| x_k - \frac{ \left\langle a_i, x_k \right\rangle}{\|a_i\|^2} a_i \|^2}
   &=  \sum_{i=1}^{m}  \frac{\|a_i\|^2}{\|A\|_F^2} \frac{ \|
   x_k - \frac{\left\langle a_i, x_k \right\rangle}{\|a_i\|^2} a_i \|^4}{\| x_k - \frac{\left\langle a_i, x_k \right\rangle }{\|a_i\|^2} a_i \|^2}\\
   &=  \sum_{i=1}^{m}  \frac{\|a_i\|^2}{\|A\|_F^2}  \left\|  x_k - \frac{\left\langle a_i, x_k \right\rangle}{\|a_i\|} \frac{a_i}{\|a_i\|} \right\|^2.
  \end{align*}
  We recall that $\|x_k\| = 1$ and thus
  \begin{align*}
   \left\|  x_k - \frac{\left\langle a_i, x_k \right\rangle}{\|a_i\|} \frac{a_i}{\|a_i\|} \right\|^2 &= \|x_k\|^2 - 2\left\langle \frac{\left\langle a_i, x_k \right\rangle}{\|a_i\|} \frac{a_i}{\|a_i\|}, x_k\right\rangle + \left\| \frac{\left\langle a_i, x_k \right\rangle}{\|a_i\|} \frac{a_i}{\|a_i\|} \right\|^2\\
   &= 1- 2\frac{\left\langle a_i, x_k \right\rangle^2}{\|a_i\|^2} + \frac{\left\langle a_i, x_k \right\rangle^2}{\|a_i\|^2} =  1- \frac{\left\langle a_i, x_k \right\rangle^2}{\|a_i\|^2}.
  \end{align*}
Altogether, we have
  \begin{align*}
  \sum_{i=1}^{m}  \frac{\|a_i\|^2}{\|A\|_F^2}  \left\|  x_k - \frac{\left\langle a_i, x_k \right\rangle}{\|a_i\|} \frac{a_i}{\|a_i\|} \right\|^2  &= \sum_{i=1}^{m}  \frac{\|a_i\|^2}{\|A\|_F^2}  \left( 1 - \frac{\left\langle a_i, x_k \right\rangle^2}{\|a_i\|^2} \right) \\
   &= \frac{1}{\|A\|_F^2} \sum_{i=1}^{m} \left(\|a_i\|^2 - \left\langle a_i, x_k \right\rangle^2\right) \\
   &= 1 - \frac{1}{\|A\|_F^2} \sum_{i=1}^{m}  \left\langle a_i, x_k \right\rangle^2 = 1 - \frac{\|Ax_k\|^2}{\|A\|_F^2}.
    \end{align*}
\end{proof}


\section{Introduction}
We study the function
\[
P_\alpha : (0,1) \to (0,1) : x \mapsto \exp(1-x^{-\alpha})
\]
with $\alpha > 0$ in the context of counting rational points. In \cite{countpointsexp}, Pila gave an upper bound on the number of rational points of bounded height on a certain exponential-algebraic surface by parametrizing it using $P_\alpha$. The bound on the derivatives of $P_\alpha$, or more generally of mild functions, are suited to bound the number of rational points on the set that it parametrizes, see \cite[Corollary 3.3]{countpointsexp}. The constant $C$ in the definition of a mild function (see Definition \ref{defmild}) plays an important role in this bound.

We will show that for the map $P_\alpha$, one can take $C = 1/\alpha$. This is sharper than the result of Pila, where he shows $C = 1+1/\alpha$ (\cite[Proposition 2.8.]{countpointsexp}) and has the advantage of being close to zero for $\alpha$ large. The main obstacle in bounding the derivatives is the polynomial in $1/x$ that appears in the computations, which is rather large.

The reason why the map $P_\alpha$ is interesting in this field, is because some functions become mild after composing with it. This has been applied successfully for $\alpha = 1$ in \cite{ccs} to construct $2$-mild parametrizations of families of subanalytic sets. In Section \ref{sec3} we show our main results. We prove such properties for some classes of functions which appear in \cite{unif} in Theorem \ref{abmmild}, from which the result on parametrizations of curves follows in Theorem \ref{curvepara}. In both results we obtain $C = 1/\alpha$ for $\alpha > 0$.

If one considers a fixed subanalytic set,  there are model theoretic results that imply that it has a $0$-mild parametrization (see for instance \cite{mildpara} for subanalytic sets). However, Yomdin has shown in \cite{yomdin2} that any analytic parametrization of the algebraic family of curves in the abstract necessarily has to contain infinitely many charts (see also \cite{ex}). More precisely, the number of charts in a $0$-mild parametrization will depend on $\epsilon$ and tend to infinity as $\epsilon$ tends to zero. We will construct an explicit $1/\alpha$-mild parametrization of this family where the size of the parametrization does not depend on $\epsilon$. Moreover, we can do this for any family of power-subanalytic curves using a pre-parametrization result of \cite{unif} or \cite{para}. Taking the above example in mind, this seems optimal. More on the field of parametrizations can be found in the survey \cite{ex}.

\section{Preliminaries}
In this short section we fix some notation for multidimensional calculus and revisit the formula for arbitrary derivatives of a function in several variables, also known as the Fa\`a di Bruno formula.

Let $f \in C^\infty(U,\R)$ be an infinitely differentiable function on an open $U \subset (0,1)^m$. For $\nu \in \N^m$ we denote $|\nu| = \nu_1+\ldots + \nu_m$ and
\[
f^{(\nu)}(x) = \left(\frac{\partial^{|\nu|}}{\partial x_1^{\nu_1}\cdots \partial x_m^{\nu_m}} f\right)(x).
\]
For $f \in C^\infty(U,\R^n)$ we set $f^{(\nu)} = (f_1^{(\nu)},\ldots,f_n^{(\nu)})$. Finally, for $x \in U$ and $\mu \in \R^m$ we denote:
\[
x^\mu = \prod_{i = 1}^m x_i^{\mu_i}.
\]

\begin{definition}[Mild functions] \label{defmild}
Let $f \in C^\infty(U,\R)$, $A,B \in \R_{>0}$ and $C \in \R_{\geq 0}$. We say that $f$ is $(A,B,C)$-mild if for any $x \in U$ and $\nu \in \N^m$:
\[
|f^{(\nu)}(x)| \leq B A^{|\nu|}(|\nu|!)^{C+1}.
\]
A map $f \in C^\infty(U,\R^n)$ is $(A,B,C)$-mild if all its component functions are. We say that $f$ is $C$-mild if it is $(A,B,C)$-mild for some $A$ and $B$.
\end{definition}

A $0$-mild function is analytic, but the converse is not true. The composition of $C$-mild functions is $C$-mild. A proof can be found in the next section. It relies on the formula for arbitrary order derivatives of a composition, which is also known as the Fa\`a di Bruno formula. In the univariate case, in propositions \ref{comp}, \ref{main} and \ref{weakpower}, we will use a more compact version stated in \cite[p.51]{dio}.

\begin{proposition}[Fa\`a di Bruno, \cite{faa}]\label{chain}
Suppose that $n$ is a positive integer, $V \subset \R^d$ and $U \subset \R^e$ are open, $f: V \to \R$, $g: U \to V$ and that $f$ and $g$ are $C^n$. For any $x \in U$ and $\nu \in \N^e$ with $|\nu| = n$ we have that: $$(f \circ g)^{(\nu)}(x) = \sum_{1 \leq |\lambda| \leq n} f^{(\lambda)}(g(x)) \sum_{s = 1}^n \sum_{p_s(\nu,\lambda)} \nu! \prod_{j = 1}^s \frac{(g^{(l_j)}(x))^{k_j}}{k_j!(l_j!)^{k_j}}$$ where $p_s(\nu,\lambda)$ is the set consisting of all $k_1,\ldots,k_s \in \N^d$ with $|k_i| > 0$ and $l_1,\ldots,l_s \in \N^e$ with $0 \prec l_1 \prec \ldots \prec l_s$ such that: $$ \sum_{i = 1}^sk_i = \lambda$$ and $$\sum_{i = 1}^s |k_i|l_i = \nu.$$ Here $l_i \prec l_{i+1}$ means that $|l_i| < |l_{i+1}|$ or, if $|l_i| = |l_{i+1}|$, then $l_i$ comes lexicographically before $l_{i+1}$.
\end{proposition}

\section{The mildness of $P_\alpha$}
In this section we prove that $P_\alpha$ is $1/\alpha$-mild. Writing $P_\alpha$ as a composition of two functions, we will use the Fa\`a di Bruno formula to bound the derivatives. The proof uses the fact that the composition of $0$-mild functions is $0$-mild. Indeed, this follows from the theory of Gevrey functions \cite{Ge}, and a proof in full generality can be found in \cite{para}. To be self contained, we briefly repeat its proof, but formulated in terms of mild functions.

\begin{proposition} \label{comp}
Suppose that $f$ and $g$ are respectively $(A_f,B_f,C)$-mild and $(A_g,B_g,C)$-mild. Then the composition $f \circ g$ is $C$-mild. Moreover, if $C = 0$ we have:
\[
A = A_g(1+A_fB_g) \quad \text{and} \quad B = \frac{A_fB_fB_g}{1+A_fB_g} < B_f.
\]
\end{proposition}

\begin{proof}
By Proposition \ref{chain}, using the triangle inequality and the mildness of $f$ and $g$, we obtain for any $x \in (0,1)$ a sum of positive terms:
\[
|(f \circ g)^{(n)}(x)| \leq \sum_{k = 1}^n B_fA_f^k (k!)^{C+1} \left( \sum_{n,k} B_n(k_1,\ldots,k_n) \prod_{i = 1}^n (B_gA_g^i(i!)^{C+1})^{k_i}\right),
\]
where the sum is over all $n$ and $k$ satisfying the relations stated in the theorem. Since $x \mapsto x^{C+1}$ is convex, taking the $(C+1)$-th rooth of the constants where needed, we are done if we can show that if $C = 0$ the sum can be bounded by $BA^nn!$.

To this end, consider the maps:
\[
\psi(x) = \frac{B_f}{1-A_f(x-B_g)} \quad \text{and} \quad \phi(x) = \frac{B_g}{1 - A_gx}.
\]
One easily verifies that $\psi^{(k)}(B_g) = B_fA_f^kk!$ and $\phi^{(k)}(0) = B_gA_g^kk!$, which are exactly the bounds on the derivatives of $f$ and $g$. Now in this case, one can compute $(\psi \circ \phi)^{(n)}(0)$ directly, without using the chain rule. One checks that:
\[
(\psi \circ \phi)^{(n)}(0) = \frac{A_fB_fB_g}{1+A_fB_g} (A_g(1+A_fB_g))^n n!.
\]
\end{proof}

Clearly, the composition of a $C_1$-mild map and a $C_2$-mild map is $\max(C_1,C_2)$-mild. We need two more elementary lemma's which will be used frequently.

\begin{lemma} \label{umild}
Suppose that $\alpha \geq 1$ and let $u_\alpha$ be the map $(0,1) \to \R$ given by $u_\alpha(x) = 1-x^{-\alpha}$. Then, for any $x \in (0,1)$ and $k > 0 \in \N$, we have:
\[
x^{\alpha+k}|u_\alpha^{(k)}(x)| \leq \alpha^k k!.
\]
\end{lemma}
\begin{proof}
This follows immediately from:
\[
u_\alpha^{(k)}(x) = (-1)^{k+1}x^{-(\alpha+k)}\alpha(\alpha+1)\cdots(\alpha+(k-1))
\]
using that $\alpha \geq 1$ to obtain the desired upper bound.
\end{proof}

Clearly, if $\alpha < 1$, one can find the upper bound $k!$ instead of $\alpha^k k!$ in the lemma.

\begin{lemma} \label{expmild}
For $r,s \in \R_{>0}$ we have:
\[
\max_{x > 0} x^{-r}\exp(-sx^{-\alpha}) = \left(\frac{r}{es\alpha}\right)^{r/\alpha}.
\]
\end{lemma}
\begin{proof}
The maximum is reached in $x = (s\alpha/r)^{1/\alpha}$.
\end{proof}

\begin{proposition} \label{main}
The map $$P_\alpha: (0,1) \to (0,1): x \mapsto \exp(1-x^{-\alpha})$$ is $(A,e,1/\alpha)$-mild, where $A = 6\alpha$ if $\alpha \geq 1$ and $A = 3 (2/\alpha)^{1/\alpha}$ for $0 < \alpha < 1$.
\end{proposition}
\begin{proof}
Defining $u_\alpha(x) = 1-x^{-\alpha}$ as in Lemma \ref{umild}, we will use the Fa\`a di Bruno formula for the map $P_\alpha = \exp \circ \, u_\alpha$. This yields:
\[
P_\alpha^{(n)}(x) = \sum_{k = 1}^n \exp(u_\alpha(x)) \left( \sum_{n,k} B_n(k_1,\ldots,k_n) \prod_{i = 1}^n (u_\alpha^{(i)}(x))^{k_i}\right).
\]
The summation over $n,k$ is subject to the conditions of Proposition \ref{chain}: $\sum_{i = 1}^n ik_i = n$ and  $\sum_{i = 1}^n k_i = k$.
Thus, in the product over the derivatives of $u_\alpha$, we obtain a factor:
\[
\prod_{i = 1}^n (x^{-(\alpha+i)})^{k_i} = x^{-\sum_{i = 1}^n (\alpha+i)k_i} = x^{-(k\alpha+n)}.
\]
Therefore, we obtain:
\begin{align*}
P_\alpha^{(n)}(x) &= \sum_{k = 1}^n x^{-(k\alpha+n)}\exp(u_\alpha(x)) \left( \sum_{n,k} B_n(k_1,\ldots,k_n) \prod_{i = 1}^n (x^{\alpha+i}u_\alpha^{(i)}(x))^{k_i}\right) \\
&=  x^{-n}e^{1-(1/2)x^{-\alpha}}\sum_{k = 1}^n x^{-k\alpha}e^{-(1/2)x^{-\alpha}} \left( \sum_{n,k} B_n(k_1,\ldots,k_n) \prod_{i = 1}^n (x^{\alpha+i}u_\alpha^{(i)}(x))^{k_i}\right).
\end{align*}
We proceed by taking absolute value and using the triangle inequality. Suppose now $\alpha \geq 1$, the other case is similar. By Lemma \ref{umild} we have that:
\[
x^{\alpha+i}|u_\alpha^{(i)}(x)| \leq \alpha^i i!
\]
and by Lemma \ref{expmild} we find that:
\[
x^{-k\alpha}e^{-(1/2)x^{-\alpha}} \leq \left(\frac{2k\alpha}{e\alpha}\right)^k \leq 2^k k!
\]
where we have used the inequality $k ^k \leq e^k k!$. Thus by Proposition \ref{comp} the sum can be bounded by $(3\alpha)^n n!$. Finally, by Lemma \ref{expmild}, we have:
\[
x^{-n}e^{1-(1/2)x^{-\alpha}} \leq e \left(\frac{2n}{e\alpha} \right)^{n/\alpha} \leq e (2/\alpha)^{n/\alpha} (n!)^{1/\alpha},
\]
using that $\alpha \geq 1$ to obtain $(2/\alpha)^{1/\alpha} \leq 2$.
\end{proof}

\section{parametrizations}
In this section we prove some results that indicate that the map $P_\alpha$ can be used to construct mild parametrizations. We start with a classical example that motivates this idea. We generalize this result in Proposition \ref{power}, which yields a parametrization result for definable families of curves. Proposition \ref{weakpower} is more general, but we obtain a $(1/\alpha + 1)$-mild map in that case.

\begin{definition}[Mild parametrization]
Let $X$ be a subset of $\R^n$. A $C$-mild parametrization of $X$ is a finite collection of $C$-mild maps $\{\phi_i: (0,1)^m \to X \mid i = 1,\ldots,N\}$ such that $X$ is covered by the images of these maps. If $X_T$ is a family of subsets of $\R^n$, then a uniform $C$-mild parametrization is a finite collection of families of $C$-mild maps $\{\phi_i: T \times (0,1)^m \to X_T \mid i = 1,\ldots,N\}$. That is: there exist $A$,$B$ and $C$ such that for any $t \in T$ the maps $\phi_{i,t}$ are $(A,B,C)$-mild and $X_t$ is covered by the union of all their images.
\end{definition}

\subsection{Yomdin's example}\label{sec41}
Consider the family of hyperbolas in $(0,1)^2$ given by $xy = \epsilon^2$, with $\epsilon \in (0,1)$. There is no uniform $0$-mild parametrization for this family (see \cite[Section 3.2]{ex}). We will show that it has a $1/\alpha$-mild parametrization consisting of three charts. We will use the map $P_\alpha$ as a power substitution, similarly to the map $x \mapsto x^2$ in \cite[Lemma 4.2]{ex} and more generally $x \to x^r$ in \cite[Proposition 4.1.5]{unif}.

First, consider the family as the graph of the map $f_\epsilon(x) = f(\epsilon,x) = \epsilon^2/x$ on the interval $\epsilon^2 < x < 1$. Now, on this interval, the first order derivative (with respect to $x$) is not uniformly bounded as $x \to \epsilon^2$. To this end we divide the domain in three parts: $\epsilon^2 < x < \epsilon$, $x = \epsilon$ and $\epsilon < x < 1$. Changing variables, the first case becomes the third case and the second case is trivial. Hence we proceed on the interval $\epsilon < x < 1$.

Let $C$ be the set $\{ (\epsilon,x) \in (0,1)^2 \mid \epsilon < x < 1\}$ and let $\tilde{C}$ be the preimage under the map $P:(0,1)^2 \to (0,1)^2: (\epsilon,x) \mapsto (P_\alpha(\epsilon),P_\alpha(x))$. Then for any $\epsilon$ and $x$ in $\tilde{C}$ we have:
\[
e^{1-\epsilon^{-\alpha}} < e^{1-x^{-\alpha}} < 1 \iff \epsilon < x < 1.
\]
We claim that for any $\epsilon$ the map $(f \circ P)_{\epsilon}$ is $1/\alpha$-mild. Moreover, the constants $A$ and $B$ do not depend on $\epsilon$ (uniform mildness). To this end, note that $$(f \circ P)(\epsilon,x) = e^{2(1-\epsilon^{-\alpha})-(1-x^{-\alpha})}.$$ We use the same approach as in the proof of the mildness of $P_\alpha$. The $n$-th derivative with respect to $x$ can be bounded in absolute value by:
\[
|(f \circ P)_{\epsilon}^{(n)}(x)| \leq \sum_{k = 1}^n e^{2(1-\epsilon^{-\alpha})-(1-x^{-\alpha})} \left( \sum_{n,k} B_n(k_1,\ldots,k_n) \prod_{i = 1}^n (|u_\alpha^{(i)}(x)|)^{k_i}\right).
\]
Writing the exponent as $2(1-\epsilon^{-\alpha}) - 2(1-x^{-\alpha}) + (1-x^{-\alpha})$ and using that $e^{(1-\epsilon^{-\alpha})} < e^{-(1-x^{-\alpha})}$, bounding the above sum completely reduces to the proof of Proposition \ref{main}. (Here we really used that we bounded the first order derivative.) Finally, one should revert the transformation in $\epsilon$ and map $(0,1)$ onto the domain of the obtained function. More precisely, denote  $\Gamma(f_\epsilon)$ for the graph of $f_\epsilon$, then we get:
\[
\phi: (0,1) \to \Gamma(f_\epsilon) : x \mapsto ( (1-P^{-1}_\alpha(\epsilon))x+P^{-1}_\alpha(\epsilon), \epsilon e^{-(1-x^{-\alpha})}).
\]
Note that we do not need the power substitution in $\epsilon$, as we will show in Theorem \ref{curvepara}.

\subsection{Results on compositions}\label{sec3}
We would like to generalize this example, using the results of \cite{unif}. Their pre-parametrization result (Theorem 4.3.1) leads to a uniform $C^r$-parametrization by using the power substitution $P_r: x \to x^r$. More precisely, given a map $\phi$ obtained by this result, the composition $\phi \circ P_r$ is mild (up to order $r$). Using the map $P_\alpha$, it might be possible to construct a uniform $1/\alpha$-mild parametrization, which would improve bounds on the number of rational points of bounded heights on sets that we can parametrize. We will show that $\phi \circ P_\alpha$ is $1/\alpha$-mild. This yields a uniform $1/\alpha$-mild parametrization for power-subanalytic curves.

We start with a more general class of functions defined in \cite{unif}, which also become mild after composing with $P_\alpha$.

\begin{definition}[Weakly mild functions]
Let $f \in C^\infty(U,\R)$, $A,B \in \R_{>0}$ and $C \in \R_{\geq 0}$. We say that $f$ is weakly $(A,B,C)$-mild if for any $x \in U$ and $\nu \in \N^m$:
\[
|f^{|\nu|}(x)| \leq \frac{B A^{|\nu|}(|\nu|!)^{C+1}}{x^\nu}.
\]
\end{definition}

The next result is similar to \cite[Proposition 4.1.5]{unif}. Intuitively, near $0$, $x \to x^r$ is a vanishing function of order $r$, thus is capable of rendering functions mild up to order $r$. In our case, this will be up to order $+\infty$. In the proof of the next proposition, we use that $\alpha \geq 1$, but clearly, it is also true for $0 < \alpha < 1$.

\begin{proposition}\label{weakpower}
Let $f:U \to \R$, $U \subset (0,1)$ open, be weakly $(A,B,0)$-mild and suppose that $f^{(1)}$ is weakly $(A,B,0)$-mild. Then the map $f \circ P_\alpha: P_\alpha^{-1}(U) \to \R$ is $(1+1/\alpha)$-mild.
\end{proposition}

\begin{proof}
Of course, we will use the Fa\`a di Bruno formula and the triangle inequality again. This leaves us with the task to bound:
\[
|(f \circ P_\alpha)^{(n)}(x)| \leq \sum_{k = 1}^n \frac{BA^kk!}{e^{(k-1)(1-x^{-\alpha})}} \left( \sum_{n,k} B_n(k_1,\ldots,k_n) \prod_{i = 1}^n (|P_\alpha^{(i)}(x)|)^{k_i}\right).
\]
Here, we have $k-1$ in the denominator because of the additional assumption on $f^{(1)}$. From the proof of Proposition \ref{main}, we know that:
\begin{align*}
|P_\alpha^{(i)}(x)| &= \left|\sum_{l = 1}^i x^{-(l\alpha+i)}e^{1-x^{-\alpha}} \left( \sum_{i,l} B_i(l_1,\ldots,l_i) \prod_{j = 1}^i (x^{\alpha+j}u_\alpha^{(j)}(x))^{l_j}\right)\right| \\
&\leq x^{-(\alpha+1)i}e^{1-x^{-\alpha}}\sum_{l = 1}^i \left( \sum_{i,l} B_i(l_1,\ldots,l_i) \prod_{j = 1}^i (x^{\alpha+j}|u_\alpha^{(j)}(x)|)^{l_j}\right).
\end{align*}
The sum over $l$ can be considered as the composition of an $(1,1,0)$-mild map and a $(\alpha,1,0)$-mild map, yielding a $(2\alpha,1,0)$-mild map. Using the relations for the sum over $n,k$ and Lemma \ref{expmild} we get:
\begin{align*}
|(f \circ P_\alpha)^{(n)}(x)| &\leq \sum_{k = 1}^n \frac{BA^kk!}{e^{(k-1)(1-x^{-\alpha})}}x^{-(\alpha+1)n}e^{k(1-x^{-\alpha})} \left( \sum_{n,k} B_n(k_1,\ldots,k_n) \prod_{i = 1}^n ((2\alpha)^ii!)^{k_i}\right) \\
&=  x^{-(\alpha+1)n}e^{1-x^{-\alpha}}\sum_{k = 1}^n BA^kk!\left( \sum_{n,k} B_n(k_1,\ldots,k_n) \prod_{i = 1}^n ((2\alpha)^ii!)^{k_i}\right) \\
&\leq e\left(\frac{(\alpha+1)n}{e\alpha}\right)^{(\alpha+1)n/\alpha} B(2\alpha(A+1))^nn! \\
&\leq eB\left(\frac{\alpha+1}{\alpha}\right)^{(\alpha+1)n/\alpha}(2\alpha(A+1))^n (n!)^{2+1/\alpha}.
\end{align*}
\end{proof}

Following the strategy of the proof of \cite[Proposition 4.1.5]{unif} for several variables, will also yield this result for weakly mild functions in more variables. Now this result is not very satisfying and it seems that it cannot be improved in general. For $\alpha = 1$ it is equal to the bound in \cite[Lemma 77]{ccs}.

However, if $f$ is of a specific form, we can improve it. The maps we consider are the maps obtained from the pre-parametrization result of \cite{unif}. Our motivating example falls within this scope. The key observation is that these maps have a good interaction with $P_\alpha$. Indeed, the main obstruction to a better result in the proof above is that we completely factor out a power of $x$ in the estimate of $|P^{(i)}_\alpha(x)|$. We will not have to do this below. Let us be more precise about which maps are given by this parametrization result.

\begin{definition}[a-b-m functions]\label{abmdef}
A map $b:U \subset (0,1)^m \to \R$ is bounded-monomial if it is either identically zero, or is of the form $x^\mu$ for some $\mu \in \R^m$ and its image is bounded. A map $b: U \to \R^n$ is bounded-monomial if all of its component functions are. A map $f: U \to \R$ is analytic-bounded-monomial (a-b-m) if $f$ is of the from
\[
f(x) = b_j(x)F(b(x))
\]
where $b$ is a bounded-monomial map, $b_j$ a component function of $b$ and $F$ is a unit on $\im(b)$, that is: it is analytic and non-vanishing on an open neighborhood of the topological closure of $\im(b)$.
\end{definition}

It is shown in \cite[Corollary 4.2.4]{unif} that if all partial derivatives of the associated bounded-monomial map $b$ of $f$ are bounded, then $f$ and all its first order derivatives are weakly $0$-mild. Thus one could apply the previous proposition to this class of functions. In the next proposition, we use that $\alpha \geq 1$ in the proof and the statement, but it is straightforward to deduce an analogue for $0 < \alpha < 1$. One just has to compute the constant $A$.

\begin{proposition} \label{power}
Suppose that $b(x) = x^{\mu}$ is a bounded-monomial map whose first order derivatives are bounded on an open $U \subset (0,1)^m$. Let $P: (0,1)^m \to (0,1)^m$ be given by
\[
P(x_1,\ldots,x_m) = (P_\alpha(x_1),\ldots,P_\alpha(x_m)).
\]
Then $b \circ P: P^{-1}(U) \to \R$ is $(A,B,1/\alpha)$-mild with $A = 2\alpha(2mN+1)$ with $N = \max_{i = 1}^m |\mu_i|$, $B = e^{|\mu|}M^2$ and $M$ a constant that depends on $b$.
\end{proposition}

\begin{proof}
We use the same strategy as in the proof of Proposition \ref{main}: using the multivariate Fa\`a di Bruno formula. We write
\[
(b \circ P)(x) = e^{\sum_{i = 1}^m\mu_i(1-x_i^{-\alpha})} = (E \circ u_\alpha)(x_1,\ldots,x_m)
\]
where $E = \exp(x_1+\ldots+x_m)$ and $u_\alpha = (u^1_\alpha(x_1),\ldots,u^m_\alpha(x_m))$ where $u^i_\alpha(x_i) = \mu_i(1-x_i^{-\alpha})$. We may suppose that all $\mu_i$ are nonzero, otherwise the map does not depend on $x_i$ and any derivative with respect to $x_i$ is trivial. We have that:
\[
(E \circ u_\alpha)^{(\nu)}(x) = \sum_{1 \leq |\lambda| \leq n} E^{(\lambda)}(u_\alpha(x)) \sum_{s = 1}^n \sum_{p_s(\nu,\lambda)} \nu! \prod_{j = 1}^s \frac{(u_\alpha^{(l_j)}(x))^{k_j}}{k_j!(l_j!)^{k_j}}.
\]
As in the proof of Proposition \ref{main}, we rewrite the formula as follows:
\[
(E \circ u_\alpha)^{(\nu)}(x) = \sum_{1 \leq |\lambda| \leq n} \frac{1}{x^{\nu+\alpha\lambda}}E^{(\lambda)}(u_\alpha(x)) \sum_{s = 1}^n \sum_{p_s(\nu,\lambda)} \nu! x^{\nu+\alpha\lambda}\prod_{j = 1}^s \frac{(u_\alpha^{(l_j)}(x))^{k_j}}{k_j!(l_j!)^{k_j}},
\]
which equals
\[
e^{|\mu|}\frac{1}{x^\nu}e^{-(1/2)\sum_{i = 1}^m \mu_ix_i^{-\alpha}}  \sum_{1 \leq |\lambda| \leq n} \frac{1}{x^{\alpha\lambda}}e^{-(1/2)\sum_{i = 1}^m \mu_ix_i^{-\alpha}}  \sum_{s = 1}^n \sum_{p_s(\nu,\lambda)} \nu! x^{\nu+\alpha\lambda}\prod_{j = 1}^s \frac{(u_\alpha^{(l_j)}(x))^{k_j}}{k_j!(l_j!)^{k_j}}.
\]
Taking the extra factor $x^{\nu+\alpha\lambda}$ into account, using the relations on $l_j$ and $k_j$, we can bound everything concerning $u_\alpha$ by the appropriate bounds of an $(\alpha,N,0)$-mild function (see Lemma \ref{umild}). Now let $I$ be such that $x_I = \min_{i = 1}^m x_i$. Let $\mu' \in \R^m$ be given by $\mu'_i = \mu_i$ for $i \neq I$ and $\mu'_I = \mu_I - 1$. Then we have:
\[
\frac{1}{x^{\alpha\lambda}}e^{-(1/2)\sum_{i = 1}^m \mu_ix_i^{-\alpha}} \leq \frac{1}{x_I^{\alpha|\lambda|}}e^{-(1/2)x_I^{-\alpha}}e^{-(1/2)\sum_{i = 1}^m \mu'_ix_i^{-\alpha}} \leq \left(\frac{2|\lambda|}{e}\right)^{|\lambda|} e^{-(1/2)\sum_{i = 1}^m \mu'_ix_i^{-\alpha}}
\]
by Lemma \ref{expmild}. Now denote $K = e^{-(1/2)\sum_{i = 1}^m \mu'_ix_i^{-\alpha}}$, then $e^{|\mu'|}K^2 = ((\partial / \partial x_I)(b)/\mu_I) \circ P$. By our assumption on the partial derivatives of $b$, taking the maximum, this can be bounded independent of $I$. Thus we find that $K < M$ for some $M$ depending on $b$ only.

Similarly we find that
\[
\frac{1}{x^\nu}e^{-(1/2)\sum_{i = 1}^m \mu_ix_i^{-\alpha}} \leq \left(\frac{2|\nu|}{e\alpha}\right)^{|\nu|/\alpha}M.
\]
Using the multivariate version of Proposition \ref{comp} (see \cite[Corollary 2.5.1]{para}) we obtain that for any $x \in U$:
\[
|(b \circ P)^{(\nu)}(x)| \leq e^{|\mu|}M^2 (2/\alpha)^{|\nu|/\alpha}(\alpha(1+m2N))^{|\nu|} (|\nu|!)^{1+1/\alpha}.
\]
\end{proof}

As we have mentioned in the example, the fact that all first order partial derivatives of $b$ are bounded is crucial. Not only for the proof, but the statement can fail if it is not bounded. Indeed, if one does not bound the first order derivative in the example of Section \ref{sec41}, one would obtain the relation $e^{2(1-\epsilon^{-\alpha})} < e^{1-x^{-\alpha}}$, which is not sufficient to bound $e^{-2(1-x^{-\alpha})}$.

\begin{theorem}\label{abmmild}
Suppose that $f: U \to \R$ is a-b-m on an open $U \subset (0,1)^m$ such that the first order derivatives of its associated bounded-monomial map $b$ are bounded. Then $f \circ P: P^{-1}(U) \to \R$ is $1/\alpha$-mild.
\end{theorem}

\begin{proof}
This follows by Proposition \ref{power} and the fact that the unit $F$ is $0$-mild since it is analytic on an open neighborhood of $\overline{\text{Im}(b)}$ \cite[Proposition 2.2.10]{primer}. This concludes the proof by the remark below Proposition \ref{comp} and the fact the the product of $C$-mild functions is $C$-mild.
\end{proof}

Before moving to our final result, let us point out that by the method above, one can show that the parametrization of the exponential-algebraic surface on the last line of \cite[p. 507]{countpointsexp} is $1/g$-mild (improving $1+1/g$). This improves the upper bound on rational points in that paper to $c(\epsilon)(\log T)^{35+\epsilon}$.

We will now use the pre-parametrization result of \cite[Theorem 3.5]{para}, which is a slight modification of \cite[Theorem 4.3.1]{unif}, to obtain a uniform $1/\alpha$-mild parametrization of a family of power-subanalytic curves. This improves the $2$-mild parametrization of subanalytic curves of \cite{ccs}.

\begin{theorem}\label{curvepara}
If $X_T$ is a power-subanalytic family of power-subanalytic curves in $(0,1)^n$, then it has a uniform $1/\alpha$-mild parametrization.
\end{theorem}

\begin{proof}
By \cite[Theorem 3.5]{para} we obtain finitely many $f: C \subset T \times (0,1) \to [-1,1]^{n-1}$, where $C$ is an open cell in $T \times (0,1)$. The map $f$ is of the following form
\[
f(x) = b_j(t,x) F(b(t,x)),
\]
where $b_j$ is a component function of $b$, all component functions of $b$ are of the form $a(t)x^r$ for some definable function $a(t)$, $r \in \R$ and $F$ is a non-vanishing analytic function on an open neighborhood of $\overline{\im(b)}$. Moreover, $(\partial/\partial x)(b)$ is (uniformly) bounded. Consider the map
\[
P: T \times (0,1) \to T \times (0,1) : (t,x) \mapsto (t,P_\alpha(x))
\]
and set $\tilde{C} = P^{-1}(C)$. It follows that for any $t$ the map $(f \circ P)_t$ is $(A,B,1/\alpha)$-mild, where $A$ and $B$ are independent of $t$ (using the same technique as in the proof of Proposition \ref{power}, using the fact that $a(t)x^r$ and its first order derivative are uniformly bounded). Since in the case of curves the walls of the cell bounding the variable $x$ only depend on $t$, for any $t$ they are automatically $(1,1,0)$-mild.
\end{proof}

It would be interesting to know whether the techniques in this paper can be applied to the cellular maps in \cite{ccs} to obtain a uniform $1/\alpha$-mild parametrization of families of subanalytic sets (of arbitrary dimension). For power-subanalytic sets, one does not really require an analogue of \cite[Proposition 4.2.6]{unif}, as mentioned in the introduction of this section. One could argue by induction instead. The obstacle in that case is the fact that $P_\alpha^{-1}(C)$ is not power-subanalytic.

\bibliographystyle{alpha}


\section*{Abstract}
Inverse probability weighting is increasingly used in causal inference, but the g-computation constitutes a promising alternative. We aimed to compare the performances of these methods with time-to-event outcomes. Given the limitations of the interpretability of the hazard ratio for causal inference, our target estimand of interest is the difference in restricted mean survival times. We report the findings of an extensive simulation study showing that both inverse probability weighting and g-computation are unbiased under a correct model specification, but g-computation is generally more efficient. We also analyse two real-world datasets to illustrate these results. Finally, we update the {\tt{R}} package {\tt{RISCA}} to facilitate the implementation of g-computation.\\

\section*{Keywords} Causal inference, Parametric g-formula, Propensity score, Restricted mean survival time, Simulation study.

\section{Introduction}
Real-world evidence is scientific evidence obtained from data collected outside the context of randomised clinical trials \cite{Sherman_RWE_2016}. The absence of randomisation complicates the estimation of the marginal causal effect (hereafter referred to merely as causal effect) of exposure (including treatment or intervention) due to a potential risk of confounding \cite{Hernan_2004}. Rosenbaum and Rubin \cite{rosenbaum_central_1983} introduced the propensity score (PS) as a tool for causal inference in the presence of measured confounders. In a binary exposure setting, it has been shown that the estimated PS is a balancing score, meaning that conditional on the estimated PS, the distribution of covariates is similar for exposed and unexposed patients. Following this property, the PS can be used in four ways to provide estimates of the causal exposure effect: matching, stratification, adjustment, and inverse probability weighting (IPW) \cite{robins_marginal_2000}. Stratification leads to residual confounding and adjustment relies on strong modelling assumptions \cite{lunceford_stratification_2004, Vansteelandt_Daniel_2014}. Although matching on PS has long been the most popular \cite{ali_reporting_2015}, IPW appears to be less biased and more precise in several studies \cite{le_borgne_comparisons_2015, Hajage__2016, austin_performance_2013}. Moreover, King and Nielsen \cite{King_PSM} argued for halting the use of PS matching for many reasons, including covariate imbalance, inefficiency, model dependence, and bias.

Causal effects can also be estimated using the g-computation (GC), a maximum likelihood substitution estimator of the g-formula \cite{Robins_1986, Snowden_2011, Keil_2014}. While IPW is based on exposure modelling, the GC relies on the prediction of the potential outcomes for each subject under each exposure status.

Several studies (see \cite{Chatton_2020} and references therein) compared the IPW and GC in different contexts. They reported a lower variance for the GC than the IPW. Nevertheless, to the best of our knowledge, no study has focused on time-to-event outcomes. The presence of right-censoring and its magnitude is of prime importance since a small number of observed events due to censoring may impact the estimation of the outcome model involved in the GC. By contrast, the IPW may perform well as long as the number of exposed patients is sufficient to estimate the PS and the total sample size is sufficiently large to limit variability in the estimated weights.
In addition to the lack of simulation-based studies comparing the GC and the IPW in a time-to-event context, the GC-based estimates, such as hazard ratio or restricted mean survival, are not straightforward to obtain. Continuous-time data result in methodological difficulties, mainly due to the time-varying distribution of individual baseline characteristics.

In the present paper, we aimed to detail the statistical framework for using the GC in time-to-event analyses. We restricted our developments to time-invariant confounders because the joint modelling of time-to-event and longitudinal data requires for controlling the time-varying confounding by GC or IPW \cite{Hernan_2001}, which is beyond the scope of this paper. We also compared the performances of the GC and IPW. The rest of this paper is structured as follows. In section 2, we detail the methods. Section 3 presents the design and findings of a simulation study. In section 4, we propose a practical comparison with two real-world applications related to treatment evaluations in multiple sclerosis and kidney transplantation. Finally, we discuss the results and provide practical recommendations to help analysts to choose the appropriate analysis method.

\section{Methods}

\subsection{Notations}

Let $(T_i,\delta_i,A_i,L_i)$ be the random variables associated with subject $i$ ($i=1,...,n)$. $n$ is the sample size, $T_i$ is the participating time, $\delta_i$ is the censoring indicator (0 if right-censoring and 1 otherwise), $A_i$ is a binary time-invariant exposure initiated at time $T=0$ (1 for exposed subjects and 0 otherwise), and $L_i=(L_{1i},...,L_{pi})$ is a vector of the $p$ measured time-invariant confounders. Let $S_a(t)$ be the survival function of group $A=a$ at time $t$, and let $\lambda_a(t)$ be the corresponding instantaneous hazard function. Suppose $D_a$ is the number of different observed times of event in group $A=a$. At time $t_j$ ($j=1,... , D_a$), the number of events $d_{ja}$ and the number of at-risk subjects $Y_{ja}$ in group $A=a$ can be defined as $d_{ja} = \sum_{i:t_i = t_j} \delta_i \mathds{1}(A_i=a)$ and $Y_{ja} = \sum_{i:t_i \geq t_j} \mathds{1}(A_i=a)$.

\subsection{Estimands}
\label{estimand}

The hazard ratio ($HR$) has become the estimand of choice in confounder-adjusted studies with time-to-event outcomes. However, it has also been contested \cite{hernan_hazards_2010, Aalen_2015}, mainly because the baseline distribution of confounders varies over time. For instance, consider a population in which there is no confounder at baseline, \textit{i.e.}, the distribution of $L$ is identical regardless of exposure level $A$. Suppose additionally that the quantitative covariate $L_1$ and the exposure $A$ are independently associated with a higher risk of death. In this situation, as illustrated in Figure \ref{L1}, the difference between the $L_1$ mean values of the two groups increases over time. Then, even when the subject-specific proportional hazard (PH) assumption is true, \textit{i.e.}, the $HR$ between an exposed and an unexposed patient with the same characteristics $L$ is constant over time, the population-average $HR$ varies over time. Instead of $HR$, one can estimate the average over time of the different time-specific $HRs$: $AHR = \int [ \lambda_1(t) / \lambda_0(t) ] f(t)$ \cite{Schemper_2009}.

Nevertheless, Aalen \textit{et al.} \cite{Aalen_2015} concluded that it is difficult to draw causal conclusions from such a relative estimand. Hern\'an\cite{hernan_hazards_2010} advocated the use of the adjusted survival curves and related differences. For instance, the restricted mean survival time (RMST) allows us to summarise a survival curve for a specific time-window and to compare two curves by looking at the difference in RMST \cite{Royston_Parmar_2013}. The RMST difference up to time $\tau$ is formally defined as :
\begin{equation}
\Delta (\tau) = \int_{0}^{\tau} [S_1(t) - S_0(t)]dt
\label{eq_rmst}
\end{equation}

This value corresponds to the difference in terms of mean event-free time between two groups of exposed and unexposed individuals followed up to time $\tau$. A further advantage of the RMST difference is its usefulness for public health decision making \cite{Poole_2010}. Note that other alternatives that might avoid this problem exist, such as the attributable fraction or the number needed to treat \cite{Sjolander_2018}.

Hereafter, we considered $AHR$ and $\Delta (\tau)$.

\subsection{Weighting on the inverse of propensity score}

Formally, the PS is $p_{i} = P(A_i=1|L_i)$, \textit{i.e.}, the probability that subject $i$ is exposed according to her/his characteristics $L_i$. In practice, analysts often use logistic regression. Let $\omega_i$ be the weight of subject $i$. Xu \textit{et al.} \cite{xu_use_2010} defined $\omega_i = A_i P(A_i=1) / p_i + (1-A_i) P(A_i=0) / (1-p_i)$ to obtain a pseudo-population in which the distribution of covariates is balanced between exposure groups, enabling estimation of the causal effect in the entire population \cite{Hernan_2004}. The use of stabilised weights has been shown to produce a suitable estimate of the variance even when there are subjects with extremely large weights \cite{robins_marginal_2000, xu_use_2010}. The weighted numbers of events and at-risk subjects at time $t_j$ in group $A=a$ are $d_{ja}^\omega = \sum_{i:t_i = t_j} \omega_i \delta_i \mathds{1}(A_i=a)$ and $Y_{ja}^\omega = \sum_{i:t_i \geq t_j} \omega_i \mathds{1}(A_i=a)$, respectively. Cole and Hern\'an \cite{cole_adjusted_2004} proposed a weighted Kaplan-Meier estimator defined as:

\begin{equation}
\hat S_a(t)= \prod_{t_j \leq t}  \Big  [ 1- d_{ja}^\omega / Y_{ja}^\omega \Big ]
\label{eq1}
\end{equation}

To estimate the corresponding $AHR$, they suggested the use of a weighted univariate Cox PH model, in which exposure is the single explanatory variable. We use equation (\ref{eq_rmst}) to estimate the corresponding $\Delta(\tau)$.

\subsection{G-computation}

Akin to the IPW, the GC involves two steps. The first consists of estimating the outcome model, called the Q-model \cite{Snowden_2011}. When suitable, it can consist of a PH regression: $h(t|A_i,L_i) = h_0(t)\exp(\gamma A_i + \beta L_i)$ where $h(t|A_i,L_i)$ is the conditional hazard function of subject $i$ at time $t$, $h_0(t)$ the baseline hazard function, and $\gamma$ and $\beta$ are regression coefficients. Estimates of the cumulative baseline hazard $\hat{H}_0(t)$ and the regression coefficients ($\hat \gamma, \hat \beta$) can be obtained by the joint likelihood approach proposed by Breslow \cite{Breslow_1972}. The second step consists of predicting the counterfactual mean survival function if all subjects would have been exposed ($do(A=1)$) or unexposed ($do(A=0)$):

\begin{equation}
\hat{S}_a(t) = n^{-1} \sum_{i=1}^n \exp \Big [- \hat{H}_0(t) \times \exp(\hat\gamma \times do(A_i=a) + \hat\beta L_i)  \Big ]
\label{eq2}
\end{equation}

Then, $\widehat{AHR}$ can be computed as the mean of the individual counterfactual hazard ratios at the observed event times \cite{Schemper_2009}:

\begin{equation}
\widehat{AHR} = \Big [ \sum_{i=1}^n \delta_i  \Big ]^{-1} \sum_{i=1}^n \delta_i \Big [\hat{\lambda}_1(t_i)/ \hat{\lambda}_0(t_i)  \Big ],
\label{eq3}
\end{equation}
where $\hat{\lambda}_a(t) = - \partial \log\hat{S}_a(t)/\partial t$, which is obtained from equation (\ref{eq2}) by numerical differentiation. We use equation (\ref{eq_rmst}) to estimate the corresponding $\Delta(\tau)$.

\subsection{Software}

We performed all the analyses using {\tt{R}} version 3.6.0. \cite{R} Source code to reproduce the results is available as Supporting Information on the journal’s web page. To facilitate their use in practice, we have implemented the previous methods in the {\tt{R}} package entitled {\tt{RISCA}} (versions $\geq$ 0.8.1), which is available at {\tt{cran.r-project.org}}. We obtained the variances by bootstrap, as recently recommended by Austin \cite{Austin_boot}.

\subsection{Identifiability conditions}

As for standard regression models, the IPW and the GC require assumptions of non-informative censoring, no measurement error, no model misspecification, and no interference \cite{Hudgens_Halloran_2008}. Three additional assumptions, called \emph{identifiability conditions}, are necessary for causal inference \cite{Naimi_2016}. (i) The values of exposure under comparisons correspond to well-defined interventions that, in turn, correspond to the versions of exposure in the data. (ii) The conditional probability of receiving every value of exposure depends only on the measured covariates. (iii) The conditional probability of receiving every value of exposure is greater than zero. These assumptions are known as \emph{consistency}, \emph{conditional exchangeability} and \emph{positivity}, respectively.

\section{Simulation study}

\subsection{Data generation}

We generated data in three steps. (a) We simulated three covariates ($L_1$ to $L_3$) from a Bernoulli distribution with parameter equal to 0.5 and three covariates ($L_4$ to $L_6$) from a standard normal distribution (Figure \ref{DAG}). (b) We generated the exposure $A$ according to a Bernoulli distribution with probability obtained by the logistic model with the following linear predictor: $-0.5 + \log(2) \cdot L_2 + \log(1.5) L_3 +\log(1.5) L_5 +\log(2) L_6$. We set the intercept to obtain the prevalence of exposed individuals at 50\%. (c) We generated the times-to-event from a Weibull PH model. We set the scale and shape parameters to 40.0 and 2.0, respectively. Based on a random variable $U_i$ drawn from a standard uniform distribution, we then computed the time-to-event from a Weibull PH model as $40.0 * [(1 - \log(1-U_i) \exp(-\gamma A_i - \log(1.8) L_1 - \log(1.3) L_2 - \log(1.8) L_4 - \log(1.3) L_5 )) - 1]^{-2.0}$, where $\gamma = \log(1.0)$ under the null hypothesis or $\log(1.3)$ under the alternative hypothesis. We subsequently censored the times-to-event using a uniform distribution on [0,70] or [0,15], leading to approximately 40\% and 90\% censored observations, respectively.

We explored three sample sizes: $n=$ 100, 500, and 2000. When the censoring rate was 90\%, we did not investigate the smallest sample size due to the reduced number of events. For each scenario, we randomly generated 10,000 datasets. To compute the difference in $\Delta (\tau)$, we defined $\tau$ in each dataset as the time at which at least 10\% of the individuals in each group (exposed or unexposed) were still at risk.

We furthermore considered two sets of covariates: $(L_1, L_2, L_4, L_5)$ the risk factors of the outcome, or $(L_2, L_5)$ the true confounders.

\subsection{Performance criteria}

We computed the theoretical values of the $AHR$ and $\Delta(\tau)$ by averaging the estimations obtained, respectively, from univariate Cox PH models ($A$ as the only explanatory covariate) and by equation (\ref{eq_rmst}) where the survival functions were estimated by the Kaplan-Meier estimator, fitted from datasets simulated as above, except $A$ was simulated independently of $L$ \cite{gayat_propensity_2012}. We reported the following criteria: (i) the percentage of datasets without convergence ; (ii) the mean absolute bias $MAB=E(\hat \theta)-\theta$, where $\theta$ is the estimand of interest; (iii) the mean square error $MSE = E[(\hat\theta - \theta)^2]$; (iv) the variance estimation bias $VEB = 100 \times (E[\widehat{SD}(\hat\theta)]/SD(\hat\theta) - 1)$, where $\widehat{SD}(\bullet)$ is the asymptotic standard deviation and $SD(\bullet)$ is the empirical standard deviation; (v) the empirical coverage rate of the nominal 95\% confidence interval (CI), defined as the percentage of 95\% CIs including $\theta$; (vi) the type I error, defined as the percentage of times the null hypothesis is rejected when the null hypothesis is true; and (vii) the statistical power, defined as the percentage of times the null hypothesis is rejected when the alternative hypothesis is true. We computed the Monte Carlo standard errors for each performance measure \cite{Morris_simulation}.

\subsection{Results}

The Monte Carlo errors were weak, and we did not encounter any convergence problems. Figures \ref{resAHR} and \ref{resRMST} present the results under the alternative hypotheses for $AHR$ and $\Delta$, respectively. The results under the null hypothesis were comparable and can be found in the supplementary material available online.

The MAB values associated with IPW and GC were similar to zero in all scenarios, regardless of the considered covariate set (all the risk factors or only the true confounders). Nevertheless, the MAB under the alternative hypothesis was slightly lower with the GC when considering all the risk factors. For instance, in terms of $\Delta$ when $n=500$ with a censoring rate of 40\%, the MAB was 0.033 for GC versus 0.061 for IPW. When considering only true confounders, these values were 0.147 and 0.066, respectively.

The GC, when considering all outcome causes, produced the best results in terms of MSE, especially for small sample sizes. For instance, when $n=100$ with a censoring rate of 40\%, the MSE related to the $AHR$ was 0.056 for GC versus 0.068 for IPW. When considering only true confounders, these values were 0.076 and 0.079, respectively.

Regarding the VEB, the results were close for the GC and IPW. One can nevertheless note lower values for the IPW when $n=100$ to estimate the $AHR$. For instance, for a censoring rate of 40\%, the VEB values were between 4.1\% and 4.7\% for GC versus 2.4\% and 3.7\% for IPW.

All scenarios broadly respected the nominal coverage value of 95\%. The power was the highest for the GC, especially when considering all the risk factors, with a gain between 10\% and 13\%.

With a censoring rate of 90\%, the two methods produced similar results, irrespective of the considered scenario.

\section{Applications}

We used data from two studies performed for multiple sclerosis and for kidney transplantation \cite{Laplaud_2019,Masset_2019}. We conducted these studies following the French law relative to clinical noninterventional research. Written informed consent was obtained. Moreover, the French commission for data protection approved the collection (CNIL decisions DR-2014-327 and 914184). To guide variable selection, we asked experts which covariates were causes of the exposure or the outcome prognosis to define the causal structure \cite{VanderWeele_Shpitser_2011}. We checked the positivity assumption and the considered covariates balance (see supplementary materials available online). The log-linearity hypothesis of continuous covariates was confirmed in the univariate analysis if the Bayesian information criterion \cite{BIC} was not reduced using natural spline transformation compared to the inclusion of the covariate in its natural scale. In case of violation, we used a natural spline transformation. We also assessed the PH assumption via the Grambsch-Therneau test at a significance level of 5\% \cite{Grambsch_Therneau_1994}. For simplicity, we performed complete case analyses.

\subsection{Dimethylfumarate versus Teriflunomide to prevent relapse in multiple sclerosis}

With the increasing number of available drugs for preventing relapses in multiple sclerosis and the lack of head-to-head randomised clinical trials, Laplaud \textit{et al.} \cite{Laplaud_2019} aimed to compare Teriflunomide (TRF) and Dimethylfumarate (DMF) using data from the multicentric cohort OFSEP. We reanalysed the primary outcome, defined as the time-to-first relapse. We presented the cohort characteristics of 1770 included patients in Table 1: 1057 patients were in the DMF group (59.7\%) versus 713 in the TRF group (40.3\%). Approximately 39\% of patients (40\% in the DMF group versus 38\% in the TRF group) had at least one relapse during follow-up.

We presented the confounder-adjusted results in the left panel of Figure \ref{forCVC}. The GC and IPW were equally robust to the considered set of covariates. For the difference in RMST, the width of the 95\% CI was lower for the GC. For instance, when we considered all the risk factors, the CI of GC had a width of 44.2 days versus 50.9 days for IPW.

The conclusion of no significant difference between TRF and DMF was unaffected by the method or the choice of the estimand.

\subsection{Basiliximab versus Thymoglobulin to prevent post-transplant complications}

Amongst non-immunised kidney transplant recipients, one can expect similar rejection risk between Thymoglobulin (ATG) and Basiliximab (BSX), two possible immunosuppressive drugs proposed as induction therapy. However, ATG may be associated with higher serious events, especially in the elderly. We aimed to reassess the difference in cardiovascular complications in ATG versus BSX patients \cite{Masset_2019}. Table 2 describes the 383 included patients from the multicentric DIVAT cohort: 204 patients were in the BSX group (53.3\%) versus 179 in the ATG group (46.7\%). Approximately 30\% of patients (29\% in the BSX group and 31\% in the ATG group) had a least one cardiovascular complication during follow-up. The median follow-up time was 1.8 years (min: 0.0; max: 8.2).

In the right panel of Figure \ref{forCVC}, we present the confounder-adjusted RMST differences for a cohort followed up to 3 years. First, the GC-based results obtained were slightly less sensitive to the considered set of covariates than those obtained by the IPW. Indeed, GC resulted in a $\hat{\Delta}(3)$ of 2.6 (IC 95\% from -1.2 to 7.0) and 4.0 (IC 95\% from 0.2 to 8.5) months with all the risk factors and only the true confounders, respectively. In contrast, the IPW resulted in a $\hat{\Delta}(3)$ of 1.5 (IC 95\% from -2.5 to 5.8) and 3.4 (IC 95\% from -0.5 to 7.9) months, respectively. Likewise, the log of $\widehat{AHR}$ obtained by GC was between -0.272 (IC 95\% from -0.703 to 0.123) and -0.440 (IC 95\% from -0.898 to -0.021) versus -0.188 (IC 95\% from -0.639 to 0.226) and -0.384 (IC 95\% from -0.830 to 0.002) for IPW. This variability in the estimations, according to the set of covariates, illustrates the importance of this choice. Second, we obtained the smallest 95\% CIs using the GC with all risk factors, in line with the corresponding lower variance previously reported in the simulations. Third, the conclusion differed depending on the method: the IPW-based results did not point to a significant difference in effect between BSX and ATG, in contrast to the GC with the true confounders.

\section{Discussion}

We aimed to explain and compare the performances of the GC and IPW for estimating causal effects in time-to-event analyses with time-invariant confounders. We focused on the average HR and the RMST difference. The results of the simulations showed that the two methods performed similarly in terms of MAB, VEB, coverage rate, and type I error rate. Nevertheless, the GC outperformed the IPW in terms of statistical power, even when the censoring rate was high. Furthermore, the simulations showed that the Q-model in the GC approach should preferentially include all the risk factors to ensure a smaller mean bias. Therefore, the main advantage of using the GC is the gain in statistical power.

In the two applications, the 95\% CIs were also narrower when using the GC with all the risk factors. Moreover, while the first application in multiple sclerosis did not highlight relevant differences between the GC and IPW, the second one in kidney transplantation illustrated the importance of the set of covariates to consider. On the basis of the GC with all risk factors, we concluded there was a significant difference between Basiliximab and Thymoglobulin. In contrast, when using the GC with only the true confounders, we reported a non-significant difference. Two arguments can explain why the consideration of all risk factors in the Q-model was associated with the best performance. First, this approach reduces the residual variance and increases the precision of the predictions of potential outcomes. Second, even if a risk factor is balanced between the exposure groups at baseline, it could become unbalanced over time (as illustrated in Figure~\ref{L1}).

Nevertheless, the higher power of the GC is counterbalanced by the inability of GC to evaluate balance of characteristics between exposure groups over time and the need for bootstrapping to estimate the variance, analytic estimators that are available for the IPW \cite{Conner_Trinquart_2019, Hajage_2018}. In practice, we must emphasise that bootstrapping the entire estimation procedure has the advantage of valid post-selection inference \cite{Efron_2014}. For instance, data-driven methods for variables selection, such as the super learner \cite{vdl_2007}, have recently been developed and may represent a promising perspective when such clinical knowledge is unavailable \cite{Blakely_2019}.

The IPW and GC are not the only available methods to estimate the causal effect. For instance, Conner \textit{et al.} \cite{Conner_Trinquart_2019} compared the performance of IPW with that of other regression-based methods. Overall, the statistical performance was similar. However, the advantage of the IPW and GC compared to other methods is the visualisation of the confounder-adjusted results in terms of the survival curve or an indicator such as RMST. The use of doubly robust estimators (DREs) could be an extension of our work \cite{Benkeser_2018, Cai_vdl_2019}. DREs combine both the PS and GC approaches for a consistent estimation of the exposure effect if a least one model has been correctly specified, circumventing the aforementioned potential model misspecification \cite{Lendle_2013}. Unfortunately, DREs can increase bias when the two models are misspecified \cite{Kang_Schafer_2007}. Moreover, several studies have reported a lower variance for the GC than DREs \cite{Chatton_2020, Kang_Schafer_2007, Lendle_2013}. The use of the GC also represents a partial solution to prevent the selection of instrumental variables \cite{Myers_2011} since it is independent of the exposure modelling.

Our study has several limitations. First, the results of the simulations and applications are not a theoretical argument for generalising to all situations. Second, we studied only logistic and Cox PH regression: other models could be applied. Keil and Edwards \cite{Keil_Edwards_2018} proposed a review of possible Q-models with a time-to-event outcome. Third, we considered only a reduced number of covariates, which could explain the abovementioned equivalence between the GC and the IPW with the extreme censoring rate. Last, we did not consider competing events or time-varying confounders that require specific estimation methods \cite{Young_2020, Daniel_2013}.

To conclude, by means simulation and two applications on real
datasets, this study tended to show the higher power of the
GC compared to IPW to estimate the causal effect with time-to-event outcomes. All the risk factors should be considered in the Q-model. Our work is a continuation of the emerging literature that questions the near-exclusivity of propensity score-based methods in causal inference.

\section*{Acknowledgements}
The authors would like to thank the members of DIVAT and OFSEP Groups for their involvement in the study, the physicians who helped recruit patients and all patients who participated in this study.  We also thank the clinical research associates who participated in the data collection. We also thank David Laplaud and Magali Giral for their clinical expertise. The analysis and interpretation of these data are the responsibility of the authors. This work was partially supported by a public grant overseen by the French National Research Agency (ANR) to create the Common Laboratory RISCA (Research in Informatic and Statistic for Cohort Analyses, {\tt{www.labcom-risca.com}}, reference: ANR-16-LCV1-0003-01) involving the development of Plug-Stat software.

\section*{Competing interest}
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Arthur Chatton obtained a grant from IDBC for this work.
Other authors received no financial support for the research, authorship, and/or publication of this article.

\section*{Supplemental material and code}
Supplemental material for this article is available online.

\bibliographystyle{vancouver}