---
abstract: |
    In many biomedical applications, outcome is measured as a “time-to-event” (eg. disease progression or death). To assess the connection between features of a patient and this outcome, it is common to assume a proportional hazards model, and fit a proportional hazards regression (or Cox regression). To fit this model, a log-concave objective function known as the “partial likelihood” is maximized. For moderate-sized datasets, an efficient Newton-Raphson algorithm that leverages the structure of the objective can be employed. However, in large datasets this approach has two issues: 1) The computational tricks that leverage structure can also lead to computational instability; 2) The objective does not naturally decouple: Thus, if the dataset does not fit in memory, the model can be very computationally expensive to fit. This additionally means that the objective is not directly amenable to stochastic gradient-based optimization methods. To overcome these issues, we propose a simple, new framing of proportional hazards regression: This results in an objective function that is amenable to stochastic gradient descent. We show that this simple modification allows us to efficiently fit survival models with very large datasets. This also facilitates training complex, eg. neural-network-based, models with survival data.

    **Keywords:** Survival Analysis, Cox Proportional Hazard model, Big data, Streaming Data, Stochastic Gradient Descent, Neural Networks.
author:
- |
    Aliasghar Tarkhan and Noah Simon\
    Department of Biostatistics\
    University of Washington\
    Seattle, WA 98195-4322, USA\
    atarkhan@uw.edu and nrsimon@uw.edu
bibliography:
- 'ref.bib'
title: 'BigSurvSGD: Big Survival Data Analysis via Stochastic Gradient Descent'
---

Introduction
============

It is commonly of interest in biomedical settings to characterize the relationship between characteristics of an individual and their risk of experiencing an event of interest [eg. progression of a disease, recovery, death, etc. see @Lee1997]. Outcomes of this type are known as “time-to-event” outcomes, and characterizing such relationships is known as Survival analysis [@Schober2018]. In such applications, we often only have partial information on some patients due to censoring (e.g., they might leave the study before experiencing the event of interest).

Cox proportional hazards regression (CoxPH) [@Cox1972] is the most common tool for conducting survival analyses. The CoxPH model assumes a particular semi-parametric relationship between the risk at each given time of experiencing an event and the features of a patient (eg. age, sex, treatment assignment, etc). To estimate, parameters in this model CoxPH regression maximizes a log-concave function known as the “partial likelihood”. Once estimated, this model can predict the person-specific risk of an event (as a function of their features). Such predictions are often used in personalized medicine, eg. in the development of prognostic and predictive biomarkers [@Bjorn2008]. To maximize the partial likelihood, it is most common to use an efficient second-order algorithm such as Newton-Raphson [@MITTAL2014] for datasets with few features (though potentially many observations). Traditionally, CoxPH has been used on data-sets with relatively few observations, though penalized extensions have been developed for high dimensional applications [@Simon2011].

It is increasingly common to have biomedical datasets with a large number of observations, especially with increasing use of electronic medical records [@Raghupathi2014]. Although the CoxPH regression model has been widely used for small-to-moderate numbers of observations, current methodologies for fitting the Cox model have issues on datasets with many observations. In particular, in fitting the Cox model, it is common to leverage the sequential structure of the partial likelihood to vastly speed up computation [from $O(n^2)$ to $O(n)$, see @Simon2011]. However, when there are a large number of observations, this can lead to computational instability (which we illustrate in this manuscript).

The second issue is that the partial likelihood does not naturally decouple over individuals or subsets of individuals. Thus, if the dataset does not fit in memory, the model can be very computationally expensive to fit: Standard distributed optimization methods such as those based the alternating direction method of multipliers [@Boyd2010] cannot be used. This additionally means that the objective is not directly amenable to stochastic gradient-based optimization methods [@Ruder2016]. Unfortunately, to fit more complex neural network-based models, it is most common to use stochastic-gradient-based optimization [@NN]. This decoupling issue makes it impossible or at least very impractical to fit neural-network-based models with time-to-event data.

In this paper, we propose a novel and simple framework for conducting survival analysis using the CoxPH model. Our framework is built upon an objective that is a modification of the usual partial likelihood function. In particular this modified objective decouples over subsets of observations, and allows us to employ stochastic-gradient-based methods that engage only a subset of our data at each iteration. We show that the parameters estimated by this new objective function are equivalent to the original parameters when the assumptions of the CoxPH model hold (and may actually be more robust in the case of model misspecification). In addition, our new objective function is amenable to optimization via stochastic-gradient based methods. Standard stochastic gradient-based algorithms are computationally efficient and stable for this objective and can easily scale to datasets that are too large to fit in memory. We discuss how our new framework can be implemented in both streaming [@Gaber2005] and non-streaming algorithms. We also discuss extending our framework to use mini-batches and we present some recommendations that we have found important, in practice, for performance.

We organize this paper as follows. In Section 2, we review the CoxPH regression model, the partial likelihood, and standard optimization tools used to maximize the partial likelihood. In Section 3, we present our new framework for fitting the CoxPH model and we prove the statistical equivalence of the parameters indicated by our optimization problem and those one from the standard CoxPH model. We also discuss applications of our proposed framework to both streaming and non-streaming algorithms. In addition, we discuss some recommendations that we have found important in practice for performance. In section \[sec4\] we discuss how our framework can be extended to left truncation and right censoring. In section \[sec5\], we discuss the equivalence of our optimization procedure to U-statistics based optimization. In Section \[Sim\], we provide simulation results that compare estimates from stochastic optimization of our modified objective to the current state of the art that estimates parameters by attempting to optimize the usual partial likelihood. In section \[Conclusion\], we conclude our paper and discuss some potential implications of our framework.

Cox Proportional Hazards (Cox PH) Model {#sec2}
=======================================

The Cox PH model proposed by Cox [@Cox1972] is a commonly used semi-parametric regression model in the medical literature for evaluating the association between the time until some event of interest and a set of variable(s) measured on a patient. More formally, suppose on each patient we measure $T$ an event time, and $X = (X_1,\ldots, X_p)$ a vector of numeric features. The Cox model engages with the so-called hazard function $$h(x,t) = \frac{p(t|x)}{S(t|x)}$$ where $p(t|x) = \frac{d}{dt}P(T < t|X=x)$, $S(t|x) = P(T > t|X=x)$. The hazard function, $h(x,t)$, can be thought of as the probability density of having an event at time $t$, given that a patient (with covariates $x$) has not had an event up until that time. In particular the Cox model assumes a particular form for the hazard function: $$\begin{aligned}
h(t, \boldsymbol{x};\boldsymbol{\beta}^*) = h_0(t)e^{f_{\boldsymbol{\beta}^*}(\boldsymbol{x})}
\label{e1}\end{aligned}$$ where $f_{\boldsymbol{\beta}^*}$ is a specified function of parameters $\boldsymbol{\beta}^*=(\beta^*_1, \beta^*_2, \dots, \beta^*_k)$ that determines the role played by $\boldsymbol{x}$ in the hazard; and $h_0(t)$ is a baseline hazard function (independent of covariates). Note that $f_{\boldsymbol{\beta}^*}(\boldsymbol{x})$ may be assumed to be of different forms in different applications: For instance, in many scenarios $k$ is taken to be $p$, and the simple linear model $f_{\boldsymbol{\beta}}(\boldsymbol{x})=\boldsymbol{x}^T\boldsymbol{\beta}=  + \beta_1 x_1 + \dots + \beta_p x_p$ is used. This model assumes that the manner in which a patient’s covariates modulate their risk of experiencing an event is independent of time. In particular it is encoded entirely in $f_{\boldsymbol{\beta}^*}$. This simplifies estimation and interpretation of the predictive model.

Our aim is to use data to estimate $\boldsymbol{\beta}^*$. In particular we will assume that we have a dataset with $n$ independent observations drawn from the model : $\mathcal{D}^{(n)}= \{\mathcal{D}_i=(y_i, \boldsymbol{x}^{(i)})|i=1,2, \dots, n\}$. For the moment we assume that there is no censoring (all event times are observed), and no ties (all event times are unique). Estimation is conducted using the log-partial-likelihood: $$\begin{aligned}
pl^{(n)}(\boldsymbol{\beta}|\mathcal{D}^{(n)}) &= \operatorname{log}\left(\prod_{i=1}^{n} \frac{h(\boldsymbol{x}^{(i)}; \boldsymbol{\beta})}{\sum_{j \in \mathcal{R}_i}h(\boldsymbol{x}^{(j)}; \boldsymbol{\beta})}\right)\nonumber\\
&= \sum_{i=1}^{n}\left(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})-\operatorname{log}\left(\sum_{j \in \mathcal{R}_i}e^{f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(j)})}\right)\right)
\label{e2}\end{aligned}$$ where $\mathcal{R}_i = \{j\,|\,t_j \geq t_i\}$ is the “risk set for patient $j$”. Note that aside from $\mathcal{R}_j$, the expression in is independent of the event times. Extending this partial likelihood to deal with censoring, left-truncation and ties is quite straightforward [@Klein2003], however for ease of exposition we do not include it in this manuscript.

Using the log-partial-likelihood from equation  an estimate of $\boldsymbol{\beta}^*$ can be obtained as $$\begin{aligned}
\hat{\boldsymbol{\beta}}^{(n)} = \underset{\boldsymbol{\beta}}{\operatorname{argmin}}\left\{-pl^{(n)}(\boldsymbol{\beta}|\mathcal{D}^{(n)})\right\}
\label{e3}\end{aligned}$$ When linear $f_{\boldsymbol{\beta}}(\bold{x})=\bold{x}^{\top}\boldsymbol{\beta}= \beta_1 x_1 + \dots + \beta_p x_p$ is used, our objective function in is convex in $\boldsymbol{\beta}$, and thus the tools of convex optimization can be applied to find $\hat{\boldsymbol{\beta}}^{(n)}$ (see *Corollary 1* in Appendix \[appA\] for the proof of convexity). In the current gold standard [survival]{} package in [R]{} [@coxph2019], Newton-Raphson is used to minimize with linear $f$. In the later sections of this manuscript we will refer to this implementation as [coxph()]{}.

For linear $f$, one can show that $\left\|\hat{\boldsymbol{\beta}}^{(n)} - \boldsymbol{\beta}^{*}\right\|_2^2 = O_p\left(n^{-1}\right)$ which is rate optimal [as is standard for estimation in parametric models, see @Van2000].

In current state-of-the-art packages, the structure of the ordered structure of the loss (as well as the gradient, and hessian) are leveraged to improve computational efficiency. In particular, we examine the gradient $$\begin{aligned}
    \nabla_{\beta}\Big\{-pl^{(n)}(\boldsymbol{\beta}|\mathcal{D}^{(n)})\Big\}= -\sum_{i=1}^{n}\Big(\dot{f}_{\boldsymbol{\beta}}(\bold{x}^{(i)})-\frac{\sum_{j \in \mathcal{R}_i}\dot{f}_{\boldsymbol{\beta}}(\bold{x}^{(j)})e^{f_{\boldsymbol{\beta}}(\bold{x}^{(j)})})}{\sum_{j \in \mathcal{R}_i}e^{f_{\boldsymbol{\beta}}(\bold{x}^{(j)})}}\Big).
    \label{e4}\end{aligned}$$ where $\dot{f}_{\boldsymbol{\beta}}(\bold{x})=\nabla_{\beta}\{f_{\boldsymbol{\beta}}(\bold{x})\}$ is the gradient of $f_{\boldsymbol{\beta}}(\bold{x})$ with respect to $\boldsymbol{\beta}$. While a naive calculate would have $n^2$ computational complexity because of the nested summations, this is not necessary. In the case that the times are ordered $t_1 < t_2 <\ldots < t_n$ we see that $R_{i} = R_{i+1}\cup \{i\}$. This allows us to use cumulative sums and differences to calculate the entire gradient in $O(n)$ computational complexity, with a single $n\operatorname{log}(n)$ complexity sort required at the beginning of the algorithm [@Simon2011]. This is also true for calculating the Hessian. Unfortunately, however, when employing this strategy, the algorithm becomes susceptible to roundoff issues, especially with a larger number of observations ($n$) and features ($p$), as seen in Section \[sec:sim-big\].

Additional inspection of the gradient in shows why stochastic-gradient-based methods cannot be used to decouple gradient calculations over observations in our sample: While the gradient can be written as a sum over indices $i=1,\ldots, n$, the denominator for the $i=1$ term involves all observations in the dataset. In the next section, we propose a novel simple modification of optimization problem  that admits an efficient stochastic-gradient-based algorithm for estimating $\beta^*$.

Big survival data analysis using SGD: BigSurvSGD {#sec3}
================================================

We begin by reformulating our problem. We consider a population parameter $\boldsymbol{\beta}^{(s)}$, defined as the population minimizer of the expected partial likelihood of $s$ random patients (which we will refer to as “strata of size s”) $$\begin{aligned}
\boldsymbol{\beta}^{(s)} = \underset{\boldsymbol{\beta}}{\operatorname{argmin}}\Big\{\mathbb{E}_{s}[-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})]\Big\}
\label{e5}\end{aligned}$$ Here we think of $\mathcal{D}^{(s)}$ as a draw of $s$ random patients from our population. Note that the minimum value for $s$ is 2, otherwise, expression (2) becomes zero for all $\beta$. By including a superscript $s$ in $\boldsymbol{\beta}^{(s)}$, we note that this parameter may depend on $s$. In fact, when the assumptions of the Cox model hold then we have $\boldsymbol{\beta}^{(s)} = \boldsymbol{\beta}^{*}$ for all $s$. The proof of this is quite simple, with details given in Appendix \[appA\].

To estimate $\beta^{*}$, we select a small fixed $s$ ($s<<n$) and directly apply stochastic gradient descent to the population optimization problem . In practice this will amount to calculating stochastic gradients using random strata of size $s$. One may note that for $s$ small, there are on the order of $n^s$ such strata. However, results for stochastic gradient descent indicate that under strong convexity of , on the order of only $n$ steps should be required to obtain a rate optimal estimator (converging at a rate of $n^{-1}$ in MSE). See Appendix \[appB\] for the proof of this convergence rate $O(n^{-1})$.

In Section \[sec:sim-big\], we see that this modification mitigates issues with roundoff error, and allows us to computationally efficiently fit survival models with millions of observations and many features.

Pairwise concordance ($s=2$)
----------------------------

An interesting special case is when we choose strata of size $s=2$, and look at pairs of patients. Then, in the case of no censoring, the population minimizer in , i.e., $\boldsymbol{\beta}^{(2)}$ maximizes the expectation of the pairwise log-partial likelihood $$\begin{aligned}
    pl^{(2)}(\boldsymbol{\beta}|\mathcal{D}^{(2)}) &= log\Big(\frac{f_{\boldsymbol{\beta}}(\bold{x}^{(1)})}{f_{\boldsymbol{\beta}}(\bold{x}^{(2)})+f_{\boldsymbol{\beta}}(\bold{x}^{(1)})}\Big)1(t_1<t_2)\nonumber\\
    &+log\Big(\frac{f_{\boldsymbol{\beta}}(\bold{x}^{(2)})}{f_{\boldsymbol{\beta}}(\bold{x}^{(1)})+f_{\boldsymbol{\beta}}(\bold{x}^{(2)})}\Big)1(t_2<t_1).
    \label{ep}\end{aligned}$$ This log-partial-likelihood can be thought of as a smoothed version of the standard concordance measure used in the concordance index [@Peter2012]. Thus, even when the proportional hazards model does not hold, the parameter $\boldsymbol{\beta}^{(2)}$ maintains a useful interpretation as the population minimizer of the average smoothed concordance index. is similar to the objective function for conditional logistic regression (CLR) with strata size $s=2$ [@Breslow1980]. In the deep learning literature, neural net models with a conditional logistic outcome layers are often referred to as Siamese Neural Networks [@Gregory2015], though, to our knowledge, these ideas have not been previously been applied to time-to-event data.

Optimization with SGD
---------------------

Suppose that we have $n_s$ independent strata, $D^{(s)}_1, \ldots, D^{(s)}_{n_s}$ each with $s$ independent patients drawn from our population (with $s\geq 2$). For ease of notation, let $I_m$ denote the indices of patients in strata $D^{(s)}_m$ for each $m \leq n_s$.

We first note, that for any $\boldsymbol{\beta}$ we have

$$\begin{aligned}
\nabla_{\boldsymbol{\beta}}\mathbb{E}_s\left[pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})\right]= \mathbb{E}_s\left[\nabla_{\boldsymbol{\beta}} \left\{pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)}_m)\right\}\right], \qquad \textrm{ for all $m\leq n_s$}
\label{e6}\end{aligned}$$

when $\bold{x}$ are drawn from a reasonable distribution (eg. bounded); and $f_{\boldsymbol{\beta}}(\bold{x})$ is not too poorly behaved (eg. lipschitz). Here $\nabla_{\boldsymbol{\beta}} \left\{pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)}_m)\right\}$ is defined analogously to using $\mathcal{D}^{(s)}_m$ $$\begin{aligned}
    \nabla_{\beta}\Big\{-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)}_m)\Big\}= -\sum_{i\in I_m} \left(\dot{f}_{\boldsymbol{\beta}}(\bold{x}^{(i)})-\frac{\sum_{j \in \mathcal{R}^m_i}\dot{f}_{\boldsymbol{\beta}}(\bold{x}^{(j)})e^{f_{\boldsymbol{\beta}}(\bold{x}^{(j)})})}{\sum_{j \in \mathcal{R}^m_i}e^{f_{\boldsymbol{\beta}}(\bold{x}^{(j)})}}\right)
    \label{e7}\end{aligned}$$ where $R_i^m = \{j\,|\,t_j \geq t_i\, \textrm{ and } i,j\in I_m\}$ are risk sets that include only patients in stratum $m$; and $\dot{f}_{\beta}$ denotes the gradient of $f$ wrt $\beta$.

From here we can give the simplest version of our stochastic gradient descent (SGD) algorithm for . We choose an initial $\hat{\boldsymbol{\beta}}(0)$ (perhaps $=0$), and at each iteration $m=1,\ldots, n_s$, we update our estimate by $$\begin{aligned}
\hat{\boldsymbol{\beta}}(m) = \hat{\boldsymbol{\beta}}(m-1) + \frac{\alpha_m}{s} \times \nabla_{\beta} \left\{pl^{(s)}(\hat{\boldsymbol{\beta}}(m-1)|\mathcal{D}^{(s)}_m)\right\}.
\label{e8}\end{aligned}$$ Here, $\frac{\alpha_m}{s}$ is the learning rate (and should be specified in advance, or determined adaptively as discussed later in Section \[sec:mmm\]). The computation time to run $n_s$ steps of stochastic gradient descent according to for linear $f_{\boldsymbol{\beta}}$ is $\sim sn_s  p = np$ (where $n = sn_s$ is the total sample size). If ordering of risk sets is not leveraged, then $\sim nps$ computation is required. In contrast, Newton’s algorithm for optimizing the full log-partial likelihood requires $\sim np^2$ computation per iteration when the roundoff-error-prone updating rule is used (and $\sim n^2p + np^2$ if not). Additionally, using these small strata of size $s$, we are not prone to roundoff issues when using stochastic optimization (because this sum is calculated separately for each strata).

### Averaging over updates

It has been shown that SGD algorithms for *strongly convex objective-functions* are asymptotically more efficient if we use a running **average** of the iterates as the final estimate [@Ruppert1988; @Polyak1992]. Additionally in this case $\alpha_m = \alpha$ can be set to a fixed value (so long as it is sufficiently small). We denote the running-average estimator by $$\begin{aligned}
\tilde{\boldsymbol{\beta}}(m) = \frac{1}{m}\sum_{i=1}^{m}\hat{\boldsymbol{\beta}}(i)
\label{e9}\end{aligned}$$ Note that the averaging process does not change the values of $\hat{\boldsymbol{\beta}}(m)$. In our simulations in the main manuscript, we use the averaged $\tilde{\boldsymbol{\beta}}(m)$. Strong convexity of the objective in depends on properties of $f_{\boldsymbol{\beta}}$ and (weakly) on the distribution of $\bold{x}$. For linear $f_{\boldsymbol{\beta}}$, and $\bold{x}$ with a non-degenerate distribution, this objective will be strongly convex (See Appendix B for the proof of strong convexity of our objective function). In such cases, standard results [@Bottou2010] show that $\left\|\tilde{\boldsymbol{\beta}}(m) - \boldsymbol{\beta}^{(s)}\right\|_2^2 = O_p(m^{-1})$. As a reminder, this is the statistically optimal rate of convergence for estimating $\boldsymbol{\beta}^{(s)}$ (or equivalently, $\boldsymbol{\beta}^{*}$ when the Cox model holds) from $m$ observations. See Appendix \[appB\] for the proof of this convergence rate $O(m^{-1})$ for averaging over iterates. When we use averaging over iterates, it has been shown that choosing the learning rate as $\alpha_m=\frac{C}{\sqrt{m}}$ (C is a constant) gives us such an optimal convergence rate [@Eric2011]. We choose this learning rate in our simulation studies in Section \[Sim\] where we tune the constant $C$ to get the desired convergence rate.

Streaming vs non-streaming
--------------------------

In the discussion above, we imagined that observations were arriving in a continual stream of strata, and that we were more concerned with the cost of computation than the cost of data collection. The algorithm we described engaged with each strata only once (in calculating a single stochastic gradient). Algorithm 1 in Appendix \[appC\] details an implementation of a streaming algorithm in this imagined scenario. In practice, we generally have a fixed (though potentially large) number of observations, $n$, that are not naturally partitioned into strata. To employ SGD here, we randomly partition our observations into $n_s$ disjoint strata of size $s$, and then carry out our updates in . Additionally, in practice we may want to take more than one pass over the data (and similarly use more than $1$ random partitioning), we discuss this further in Section \[sec:mmm\]. Algorithm 2 in Appendix \[appC\] details an implementation that takes multiple passes over the data (and includes additional bells and whistles discussed in Section  \[sec:mmm\]).

### Mini-batches, Moment-based-learning-rate, and Multiple Epochs {#sec:mmm}

For small $s$, the stochastic gradients given in are potentially quite noisy. In these cases, rather than using only a single stratum to calculate the stochastic gradient, it may be preferable to average results over multiple strata. This is known as a using a mini-batch [@Ruder2016] in the SGD literature. In this survival context, we use batches of strata, so a larger strata size $s$ can eliminate the need for mini-batch sizes of greater than $1$.

Choosing a reasonable value for the hyper-parameter $\alpha_m$ (the learning rate) is critical for good practical performance of the algorithm. A number of publications have developed methods for adaptively selecting the learning rate using moments, including [Adagrad]{} [@Duchi2011], and [ADAM]{} [@Kingma2014]. However, it was shown that [ADAM]{} may not converge in some settings and an updated version called [AMSGrad]{} has been proposed [@Sashank2019]. We found substantially improved performance on simulated data with [AMSGrad]{} over the simple non-adaptive updating rule given in , especially in combination with averaging as discussed above.

In practice, taking only a single pass over the data leads to poor empirical performance. We generally use multiple passes (or “epochs”). In particular, for each epoch, we (1) randomly partition our data into strata; then (2) use the last updated iterate of $\hat{\boldsymbol{\beta}}$ from the previous epoch as the initial iterate of $\hat{\boldsymbol{\beta}}$ for the new epoch; and finally (3) apply a full pass of stochastic gradient descent (with eg. averaging and momentum) over the partitioned data. In practice we found that $\sim 100$ epochs was more than enough for very robust convergence (often $10$ was fine).

Extensions for Left Truncation and Right Censoring {#sec4}
==================================================

In practice it is common for participants in a study to leave the study before experiencing an event. This is known as right censoring. In particular it is often assumed that a patient has some random “time until censoring” $C$, and what we observe is $\operatorname{min}(C,T)$, whichever happens first (the event or censoring), along with $\delta = I(T<C)$ an indicator that the patient experienced an event (rather than censoring). Censored patients still contribute some information to estimation of $\boldsymbol{\beta}^{*}$ — in particular, if a patient is censored quite late, then there is a long period of time during which we know they did not experience an event (so likely we should estimate them to be low risk).

In some studies it is common to consider the event time $T$ as the age at which a patient had an event (rather than the calendar on study). This means that patients do not enroll in the study at $T=0$ (and patients can only be observed to fail once they are enrolled). This phenomenon, wherein patients enroll in a study at times other than $T=0$, is known as left-truncation.

When i) the assumptions of the Cox model hold; and ii) censoring and truncation times are independent of event times conditional on covariates $\bold{x}$, it is relatively straightforward to adapt cox regression to accommodate these missingness mechanisms. The partial likelihood is modified in 2 minor ways: 1) The outer summation is only taken over indices for which an event occurred (i.e. that were not censored); and 2) The risk sets, $R_i$, are modified to include only patients currently at risk (who have already been enrolled, and have not yet been censored, or had an event) [@Andersen1993]. We can similarly modify our objective function , and apply our algorithm with only minor modification (a slight change in the gradient).

U-statistic based optimization {#sec5}
==============================

While in this manuscript we discuss obtaining an estimator by directly attempting to minimize the population objective function using SGD, there is a corresponding empirical minimization problem. In particular, for strata of size $s$, one could define an estimator $\hat{\beta}^s$ by $$\hat{\beta}^{(s)} = \underset{\boldsymbol{\beta}}{\operatorname{argmin}}\Big\{-\sum_{\mathcal{D}^{(s)} \subset \mathcal{D}^{(n)}}pl^s(\boldsymbol{\beta}|\mathcal{D}^{(s)})\Big\}
\label{e10}$$ As in a standard U-statistic [@Van2000], this sum is taken over all subsets of $s$ patients out of our original $n$ patients (resulting in ${n\choose s} \sim n^s$ terms). This appears to be a difficult optimization problem, given the enormous number of terms. However, our approach shows that, in fact, only $\sim n$ of those terms need be considered: The majority contain redundant information. In fact, one could see this directly by noting that the objective function in  decouples over subsets: An application of stochastic gradient descent here would involve sampling strata with replacement, and (assuming the objective is strongly convex), with averaging over iterates would converge to a tolerance of $1/n$ after $n$ steps. This approach which is known as *incomplete U-statistics* [@Gunnar1976] could be taken more generally for losses defined by $U$-statistics.

Simulation Results {#Sim}
==================

**Data Generation**

We assume that our event time follows the Cox PH model with simple linear $f_{\boldsymbol{\beta}}$ detailed below. We generate the baseline hazard $h_0(t)$ using an exponential distribution with parameters $\lambda=1$. We generate the censoring and event times independently. The details of the data simulation procedure are given below [@Bender2005] $$\begin{aligned}
&\boldsymbol{\beta^*}=\boldsymbol{1}_{P\times 1}\\
&X_i \sim Uniform(-\sqrt{3},\sqrt{3})\,\, (\text{unit-variance uniform variable}),\\
&y_i \sim exp(\mu=e^{-\boldsymbol{X}_i^T\boldsymbol{\beta^{*}}})\,\, (\text{time to event/censoring}),\\
&\delta_i \sim Bernoulli(p=1-p_c),\,\, p_c = Pr(t_i>c_i)\\
\label{e12}\end{aligned}$$ where $y_i=min(t_i, c_i)$, i.e., time to event or censoring whichever comes first. Here $p_c$, the probability of censoring, is a parameter we can tune. Though this is not written in the form of , it is still consistent with the Cox PH model assumptions, with $f_{\boldsymbol{\beta^{*}}}(\boldsymbol{x}) = {\boldsymbol{\beta^{*}}}^{\top}\boldsymbol{x}$. In all comparisons, we include the performance of [coxph()]{} the gold standard [R]{} implementation of Newton’s algorithm for maximizing the partial likelihood.\
We used R version 3.6.1 [@R2019] to conduct the analyses. All simulations were conducted on a quad-core Intel Core i7-3520 M CPU @ 2.9GHz with 12 GB RAM. In all figures, we used mean of mean-square-error (MSE) over datasets for central tendency (curves) and standard error of MSE over datasets for variability (error bars o curves). Our implementation of the SGD procedure described in this manuscript is publicly available in the github repository: [https://github.com/atarkhan/bigSurvSGD]{} [@Tarkhan2020]. The implementation is written in [R]{}, with the computational back-end written in [C++]{}.

Small Data Results {#small-data-results .unnumbered}
------------------

We first evaluate the statistical efficiency of our estimation procedure (using strata sizes of less than $n$). We evaluate mean-squared error (MSE) between estimators $\tilde{\beta}$, and $\hat{\beta}^{(n)}$ and the truth, $\beta^{*}$ over 1000 simulated datasets with number of epochs up to 100. We see that for small strata sizes (eg. 2 and 5), there is some statistical inefficiency: Even though the convergence rate is still $\frac{1}{n}$, the constant in front appears to be worse with small strata sizes. For large strata sizes, there appears to be nearly no statistical inefficiency. For practical purposes, the SGD-based estimators do quite well. This can be seen in Figure \[fig01\].

We next evaluate the performance of averaged SGD with a fixed learning rate, against averaged SGD with an adaptive learning rate (using [AMSGrad]{}) with a fixed strata-size of $20$ for both. In addition, we try various numbers of epochs (from 10 to 100). Performance is shown in Figure \[fig02\] for 1000 simulated datasets. We see that with enough epochs (around 100) both perform well. However, [AMSGrad]{} nearly reaches that performance with as few as 50 epochs, where using a fixed learning rate does not attain that performance with fewer than $100$ epochs. For both of these methods, we tuned our \[initial\] learning-rate to be optimal in these experiments.

![$\operatorname{log}_{10}$(MSE) of estimates from AveAMSGrad (with optimal $C$ for the learning rate defined as $\alpha_m=\frac{C}{\sqrt{m}}$) for different strata sizes ($S$) with B=1, P=10, 100 epochs, and probability of censoring $p_c=0.2$.[]{data-label="fig01"}](fig1BarNoTitle.pdf "fig:")\

![$\operatorname{log}_{10}(MSE)$ of estimates from averaged SGD (AveSGD), averaged AMSGrad (AveAMSGrad), and . We choose the proportionality constant $C=3$ and $C=2$ for learning rate defined as $\alpha_m=\frac{C}{\sqrt{m}}$ for AveAMSGrad and AveSGD, respectively), S=20, B=1, P=10, and probability of censoring $p_c=0.2$.[]{data-label="fig02"}](fig2BarNoTitle.pdf)

In practice, we have found that [AMSGrad]{} is much more robust to misspecification of this initial learning rate. Figure \[fig03\] compares MSE of the estimate from AveAMSGrad for different choices of the proportionality constant in the learning rate over 1000 simulated datasets. We see that selecting the constant $C$ around 2$\sim$5 in our learning rate (defined as $\alpha_m=\frac{C}{\sqrt{m}}$) gives strong performance. However, AveAMSGrad is relatively robust to a wide range of the learning rates around the optimum value due to its capability of adapting the learning rate over iterates.

![$\operatorname{log}_{10}$(MSE) of estimates from AveAMSGrad based on sample size ($n$) with $B=1$, $P=10$, $S=20$, probability of censoring $p_c=0.2$, and different proportional constant $C$ for learning rate defined as $\alpha_m=\frac{C}{\sqrt{m}}$[]{data-label="fig03"}](fig3BarNoTitle.pdf)

We discussed the computational instability of [coxph()]{} in small-to-moderate sized datasets and how our framework can avoid such an instability. Figure \[fig04\] compares the MSE of estimates from for the small-to-moderate sample size ($n$) and number of features ($P$) over 1000 simulated datasets. As we see, performs poorly for larger $P$ and $n$. For instance, with $(P=50, n=1000)$ performs worse than $(P=50, n=100)$. This is because of computational instability with for the larger sample sizes and numbers of features. One important feature of these examples is that we include a large amount of signal (which increases as the number of features increases). With less signal, this instability is less pronounced unless very large sample sizes are used.

![$\operatorname{log}_{10}(MSE)$ of estimates from for different sample size ($n$) and different dimension of coefficient ($P$) with $s=20$ and probability of censoring $p_c=0.2$.[]{data-label="fig04"}](Fig4BarNoTitle.pdf)

Big Data Results {#sec:sim-big .unnumbered}
----------------

We next consider the numerical stability of our algorithm/framework (versus directly maximizing the full partial likelihood using Newton’s algorithm). We generated 100 datasets and we only used one epoch for AveAMSGrad algorithm. Figure \[fig05\] shows a surprising and unfortunate result for [coxph()]{}: We see that as sample size increases drastically, the performance of [coxph()]{} actually stops improving and starts getting worse! In particular for $p = 20$, [coxph()]{} is basically producing nonsense by the time we get to $10,000$ observations for this simulation setup. This indicates that for large datasets the current gold standard may be inadequate, though we do note that there is a large amount of signal in these simulations (more than we might often see in practice). In contrast [BigSurvSGD]{} has no such issues and gives quite strong performance even with only one epoch. Note did not add error bars for ease of illustration for this case.

![MSE of estimates from AveAMSGrad (with proportional constant $C=3$ for learning rate defined as $\alpha_m=\frac{C}{\sqrt{m}}$) and for different values of sample size ($n$) and dimension ($P$) with B=1, s=20, and probability of censoring $p_c=0.2$.[]{data-label="fig05"}](fig5BarNoTitle.pdf "fig:")\

We next examine the computational efficiency of our framework and algorithm. Figure \[fig06\] plots computing time (in seconds) for estimates of from and AveAMSGrad for different sample sizes (we did not add error bars for ease of illustration). We considered four different numbers of covariates $P=10, 20, 50$ and $P=100$ to examine the sensitivity of the computing time to the dimension of $\boldsymbol{\beta}$. In these examples, our algorithm read the data in chunks from the hard-drive (allowing us to engage with datasets difficult to fit in memory). The computing time of our proposed framework increases linearly in the sample size, $n$, (and in $p$). The computing time for also grows roughly linearly in $n$, though it grows quadratically in $p$ (which can be seen from the poor performance with $p=100$). Furthermore, fails for the medium-to-large datasets as it is poorly equipped to deal with datasets that do not easily fit in memory ([R]{} unfortunately generally deals somewhat poorly with memory management). For instance from Figure \[fig06\], with $P=100$ crashes after $10^5$ observations. As a reminder, the statistical performance of the output of [coxph()]{}, due to floating-point issues, degenerates much earlier.

![Computation time (seconds) of estimates from [coxph]{} and AveAMSGrad (with C=3 for $\alpha_m=\frac{C}{\sqrt{m}}$) based on sample size (n) with B=1, s=20, probability of censoring $p_c=0.2$.[]{data-label="fig06"}](Fig2bNoTitle.pdf "fig:")\

Conclusion {#Conclusion}
==========

We propose a simple and novel framework for doing survival analysis under a Cox proportional hazards model. Our framework leverages a modified optimization problem which allows us apply iterative methods over only a subset of our observations at each time. In particular it allows us to leverage the tools of stochastic gradient descent (and its extensions). This results in an algorithm that is more computationally efficient and stable than the current state of the art. We showed that our framework can handle large survival datasets with little difficulty. This framework will also facilitate training complex models such as neural network with large survival datasets.

Appendix {#appendix .unnumbered}
========

Proof of consistency for parameter $\bold{\beta}^{(s)}$ in section 3 {#appA}
====================================================================

\[app:theorem\]

In this appendix, for the sake of completeness, we prove the Fisher consistency of parameter $\boldsymbol{\beta}^{(s)}$. We treat the Cox proportional hazard model as a counting process [@SurvPoint2008]. We assume that censoring and survival times are independent given the covariate vector of interest $\boldsymbol{x}$ and they follow the model (1) with true parameter $\boldsymbol{\beta}^*$. In the following, we define some terminology before proceeding with the proof.

***Definition 1: $dN_i(u)$***. For patient $i$ with time to event $t_i$ define the counting process $dN_i(u)$ by $$\begin{aligned}
    \int_a^b g(u)dN_i(u) = \Bigg\{ {0\quad\quad\,\,\, if\quad t_i\not\in [a,b] \atop g(t_i)\quad if\quad t_i\in [a,b]}.\end{aligned}$$ For instance, if we define $g(u)=1$, the above expression is an indicator representing whether patient $i$ failed in interval $[a,b]$ (i.e., 1 represents failure and 0 otherwise). We further define $dN^{(s)}(u)=\sum_{i}^s dN_i(u)$ which is a counting process for failure times over all $s$ patients. We assume that the failure time process is absolutely continuous w.r.t. Lebesgue measure on time so that there is at most one failure at any time $u$ (i.e., no ties).

***Definition 2: $M_i(u)$***. We define $M_i(u)$ to be an indicator representing whether patient $i$ is at risk at time $u$, i.e., $t_i \leq u$. By this definition, $M(u)=\sum_{i}^sM_i(u)$ indicates number of patients who are at risk at time $u$. Note that the independent censoring assumption implies that those $M(u)$ patients at risk at time $u$ (who have not yet failed or been censored) represent a random sample of the sub-spopulation of patients who will survive until time $u$.

***Definition 3: $F(u)$***. Let $F(u)$ denote the filtration that includes all information up to time $u$, i.e., $$\begin{aligned}
    F(u)=\{(dN_i(t), M_i(t), x^{(i)}), i=1, \dots, s\,\quad \text{for $t<u$ and $dN^{(s)}(u)$}\}\end{aligned}$$ Note that given $F(u)$, we know whether patient $i$ failed or was censored (i.e., $\delta_i$), when they failed/were censored (i.e., $y_i$), and their covariate vectors $\boldsymbol{x}^{(i)}$.\
Using the above definitions, one can write the log-partial-likelihood as $$\begin{aligned}
    pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)}) = \sum_{i=1}^s \int_0^\tau \Big\{f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})-log\Big(\sum_{l=1}^s M_l(u)exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}))\Big)\Big\} dN_i(u)\end{aligned}$$ where $\tau$ is the duration of the study. Note that we keep $f_{\boldsymbol{\beta}}(\boldsymbol{x})$ in the most general form and we only assume that it is differentiable in $\boldsymbol{\beta}$ for all $\boldsymbol{x}$. Then the score function may be written as $$\begin{aligned}
    U^{s}(\boldsymbol{\beta}) &= \frac{\partial pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})}{\partial \boldsymbol{\beta}} \nonumber\\
    &= \sum_{i=1}^s \int_0^\tau \Big(\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})-\frac{\sum_{l=1}^s M_l(u)\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)})exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}))}{\sum_{l=1}^s M_l(u)exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}))} \Big)dN_i(u)\nonumber\\
    &= \sum_{i=1}^s \int_0^\tau \Big(\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})-\sum_{l=1}^{s}w_l\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}) \Big)dN_i(u)\nonumber\\
    &= \sum_{i=1}^s\int_0^\tau \dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})dN_i(u)-\sum_{l=1}^{s}\int_0^\tau w_l\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}) dN^{(s)}(u).
    \label{A2}\end{aligned}$$ where $w_l = \frac{M_l(u)exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}))}{\sum_{l=1}^s M_l(u)exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}))}$ is a weight proportional to the hazard of failure of patient $l$; $dN^{(s)}(u)=\sum_{i=1}^sdN_i(u)$.\
Now we show that the parameter $\boldsymbol{\beta^{(s)}}$ is Fisher consistent, i.e., $E[U^{s}(\boldsymbol{\beta^*})]=0$. $$\begin{aligned}
    E_{\boldsymbol{\beta^*}}[U^{s}(\boldsymbol{\beta^*})] &= E_{\boldsymbol{\beta^*}}\Big[\sum_{i=1}^s\int_0^\tau \dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})dN_i(u)-\sum_{l=1}^{s}\int_0^\tau w_l\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)})dN^{(s)}(u)\Big]\nonumber\\
    &=\sum_{i=1}^s\int_0^\tau E_{\boldsymbol{\beta^*}}\Big[\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})dN_i(u)\Big]-\sum_{l=1}^{s}\int_0^\tau E_{\boldsymbol{\beta^*}}\Big[w_l\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)})dN^{(s)}(u)\Big]\nonumber\\
    &\stackrel{(a)}{=}\sum_{i=1}^s\int_0^\tau E_{F(u)}\Big[E_{\boldsymbol{\beta^*}|F(u)}[\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})dN_i(u)]\Big]\nonumber \\
    &\,\,- \sum_{l=1}^{s}\int_0^\tau E_{F(u)}\Big[E_{\boldsymbol{\beta^*}|F(u)}[w_l\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)})]dN^{(s)}(u)\Big]\nonumber\\
    &\stackrel{(b)}{=}\sum_{i=1}^s\int_0^\tau E_{F(u)}\Big[E_{\boldsymbol{\beta^*}|F(u)}[\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})dN_i(u)]\Big]\nonumber\\ 
    &- \sum_{l=1}^{s}\int_0^\tau E_{F(u)}\Big[w_l\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)})dN^{(s)}(u)\Big]\nonumber\\
    &\stackrel{(c)}{=}\sum_{i=1}^s\int_0^\tau E_{F(u)}\Big[w_i\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})dN^{(s)}(u)\Big]\nonumber\\ 
    &- \sum_{l=1}^{s}\int_0^\tau E_{F(u)}\Big[w_l\dot{f}_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)})dN^{(s)}(u)\Big]\nonumber\\
    &= 0
    \label{A3}\end{aligned}$$ where (a) follows from the conditional expectation given the filtration $F(u)$; (b) follows from the fact that $w_l$ and $\dot{f}_{\boldsymbol{\beta^*}}(\boldsymbol{x}^{(l)})$ are known given $F(u)$; (c) follows from the fact that $E_{\boldsymbol{\beta^*}|F(u)}[dN_i(u))]=w_idN(u)$. Since $E_{\boldsymbol{\beta^*}}[U^{s}(\boldsymbol{\beta^*})]=0$, the parameter $\boldsymbol{\beta}^{(s)}$ is Fisher consistent. In the following, we present a sufficient condition under which such a parameter is a global minmizer of $E[-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})]$.\
***Corollary 1:*** The $\boldsymbol{\beta^{(s)}}$ is a global minimizer of $E[-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})]$ if $f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})$ is an affine function of $\boldsymbol{\beta}$.\
*Proof*. Suppose that we choose $f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)})$ a convex function of $\boldsymbol{\beta}$. Then, the first term inside summation of $-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})$ (i.e., $-f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)})$) is a concave function. We show that the second term inside summation of $-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})$, i.e., $log\Big(\sum_{l=1}^s Y_l(u)exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}))\Big)$ is a convex function through the following steps:\
*Step 1:* $exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}))$ is a convex function because $exp(.)$ is a non-decreasing convex function and $f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)})$ is convex [@Boyd2004 sec. 3.2.4].\
*Step 2:* $\sum_{l=1}^s Y_l(u)exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}))$ is a convex function because non-negative weighted sums of convex functions is convex [@Boyd2004 sec. 3.2.1].\
*Step 3:* $log\Big(\sum_{l=1}^s Y_l(u)exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}))\Big)$ is a convex function because $log(.)$ is non-decreasing concave function and $\sum_{l=1}^s Y_l(u)exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)})$ is a convex function [@Boyd2004 sec. 3.2.4].\
Therefore, expression $-f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})Y_i(u)+log\Big(\sum_{l=1}^s Y_l(u)exp(f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(l)}))\Big)$ is the sum of a concave function and a convex function. A sufficient condition for convexity of this expression is convexity of $-f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})$ or equivalently concavity of $f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})$. This means that $f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})$ needs to be both convex and concave at the same time. The only functions that satisfy this condition are affine functions [@Boyd2004 sec. 3.1.1]. By choosing $f_{\boldsymbol{\beta}}(\boldsymbol{x}^{(i)})$ as an affine function, $-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})$ and hence $E[-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})]$ become convex functions. Therefore, the parameter $\boldsymbol{\beta^{(s)}}$ becomes a global minimizer of $E[-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})]$. $\square$

Having the loss function $E[-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})]$ a convex function motivates us to explore the convergence rate of SGD-based minimization algorithms in the next section.

Convergence rate of SGD-based estimate {#appB}
======================================

In this appendix, for the sake of completeness, we prove that our SGD-based estimate can achieve the convergence rate of $O(n^{-1})$ for choices $f_{\boldsymbol{\beta}}(\bold{x})=\boldsymbol{\beta}^T\bold{x}$, $s=2$, and no ties. We also assume that our covariates are bounded, i.e., there exists some $C < \infty$ such that $||\bold{x}||<C$ with probability $1$ and that we consider a domain of optimization $\bold{B}$ such that $\underset{\boldsymbol{\beta}\in \bold{B}}{\operatorname{max}}\Big\{||\boldsymbol{\beta}^T\boldsymbol{x}||\Big\}<\infty$.\
For the sake of simplicity of proof, we rewrite our optimization procedure as $$\begin{aligned}
\boldsymbol{\beta}^{(s)} &= \underset{\boldsymbol{\beta}}{\operatorname{argmin}}\Big\{\mathbb{E}_{s}[-pl^{(s)}(\boldsymbol{\beta}|\mathcal{D}^{(s)})]\Big\}\nonumber\\
&= \underset{\boldsymbol{\beta}}{\operatorname{argmin}}\Big\{\mathbb{E}_{s}[L^{(s)}(\boldsymbol{\beta})]\Big\}\nonumber\\
&= \underset{\boldsymbol{\beta}}{\operatorname{argmin}}\Big\{L(\boldsymbol{\beta})\Big\},
\label{B1}\end{aligned}$$ and we estimate $\boldsymbol{\beta}^{(s)}$ iteratively using SGD from strata of size $s$ through $$\begin{aligned}
\hat{\boldsymbol{\beta}}(m) = \hat{\boldsymbol{\beta}}(m-1) - \alpha_m \times\nabla_{\boldsymbol{\beta}}L^{(s)}(\hat{\boldsymbol{\beta}}(m-1)|\mathcal{D}^{(s)}_m).
\label{B2}\end{aligned}$$ Authors in [@Eric2011] showed that if the loss function $L^{(s)}(\boldsymbol{\beta})$ satisfies the 4 following conditions, then, we can achieve the optimal convergence rate of $O(m^{-1})$ by choosing $\alpha_k=Cm^{-1}$ for the single SGD-based iterates $\hat{\boldsymbol{\beta}}(m)$ and $\alpha_k=Cm^{-1}$, $r \in [0.5, 1)$ for averaging over iterates $\tilde{\boldsymbol{\beta}}(m)$ (i.e., Polyak-Ruppert average).\
**Condition 1:** Gradient of $L^{(s)}(\boldsymbol{\beta})$ is an unbiased estimate of the gradient $L(\boldsymbol{\beta})$,\
**Condition 2:** Gradient of $L^{(s)}(\boldsymbol{\beta})$ is D-Lipschitz-continuous, i.e., $\forall$ $k\geq 1$ and $\forall$ $\boldsymbol{\beta}_1, \boldsymbol{\beta}_2 \in \mathbb{R}^p$, there exists $D\geq0$ such that, $$\begin{aligned}
    ||\nabla_{\boldsymbol{\beta}}L^{(s)}(\boldsymbol{\beta}_1)-\nabla_{\boldsymbol{\beta}}L^{(s)}(\boldsymbol{\beta}_2)||\leq D||\boldsymbol{\beta}_1-\boldsymbol{\beta}_2||,
    \label{B3}\end{aligned}$$

**Condition 3:** $L(\boldsymbol{\beta})=\mathbb{E}_{s}[L^{(s)}(\boldsymbol{\beta})]$ is $\mu$-strongly convex, i.e., $\forall$ $\boldsymbol{\beta}_1, \boldsymbol{\beta}_2 \in \mathbb{R}^p$, there exists $\mu>0$ such that, $$\begin{aligned}
    L(\boldsymbol{\beta}_1) \geq L(\boldsymbol{\beta}_2)+\nabla_{\boldsymbol{\beta}}L(\boldsymbol{\beta}_2)^T(\boldsymbol{\beta}_1-\boldsymbol{\beta}_2)+\frac{\mu}{2}||\boldsymbol{\beta}_1-\boldsymbol{\beta}_2||^2,
    \label{B4}\end{aligned}$$

**Condition 4:** Variance of the gradient of $L^{(s)}(\boldsymbol{\beta})$ is bounded, i.e., there exists $\sigma^2 \in \mathbb{R}_{+}$ such that for all $k\geq$, $$\begin{aligned}
    \mathbb{E}(||\nabla_{\boldsymbol{\beta}}L^{(s)}(\boldsymbol{\beta}^*)||^2) \leq \sigma^2, w.p.1,
    \label{B5}\end{aligned}$$

In the following, we show that the loss function in our framework, i.e., $L^{(s)}(\boldsymbol{\beta})$ satisfies all four conditions above.\
**Proof of Condition 1:**\
This condition is automatically satisfied based on the definition of $L(\boldsymbol{\beta})=\mathbb{E}_{s}[L^{(s)}(\boldsymbol{\beta})]$ in and that $\nabla_{\boldsymbol{\beta}}L(\boldsymbol{\beta})=\nabla_{\boldsymbol{\beta}}\mathbb{E}_{s}[L^{(s)}(\boldsymbol{\beta})]=\mathbb{E}_{s}[\nabla_{\boldsymbol{\beta}}L^{(s)}(\boldsymbol{\beta})]$.\
**Proof of condition 2:** The loss function $L^{(s)}(\boldsymbol{\beta})$ belongs to $C^\infty$ continuous function family and proving is equivalent to proving $$\begin{aligned}
    \exists\, D\geq0,\, s.t., \,\forall \nu \in S_\nu=\{\nu: ||\nu||_2 = 1\},\quad \nu^T\nabla^2_{\boldsymbol{\beta}}L^{(s)}(\boldsymbol{\beta})\nu \leq D
    \label{B6}.\end{aligned}$$ For $f_{\boldsymbol{\beta}}(\boldsymbol{x})=\boldsymbol{\beta}^T\boldsymbol{x}$, $s=2$ and assuming no ties, $\nabla^2_{\boldsymbol{\beta}}L^{(s)}(\boldsymbol{\beta})$ can be simplified as $$\begin{aligned}
    \nabla^2_{\beta}\Big\{L^{(s)}(\boldsymbol{\beta})\Big\}&= w(1-w)\bold{X}\begin{bmatrix}
1 & -1 \\
-1 & 1
\end{bmatrix}\bold{X}^T\nonumber\\
&\times\Big(1(\delta_1=1)1(y_1<y_2)+1(\delta_2=1)1(y_2<y_1)\Big)\nonumber\\
&=w(1-w)(\bold{x}^{(1)}-\bold{x}^{(2)})(\bold{x}^{(1)}-\bold{x}^{(2)})^T\nonumber\\
&\times\Big(1(\delta_1=1)1(y_1<y_2)+1(\delta_2=1)1(y_2<y_1)\Big)
    \label{B7}\end{aligned}$$ where $\bold{X} = [\bold{x}^{(1)}, \bold{x}^{(2)}]$; $w = w_1(\bold{x}^{(1)},\bold{x}^{(2)}, \boldsymbol{\beta})= 1-w_2(\bold{x}^{(1)},\bold{x}^{(2)}, \boldsymbol{\beta})=\frac{e^{\boldsymbol{\beta}^T\boldsymbol{x}^{(1)}}}{e^{\boldsymbol{\beta}^T\boldsymbol{x}^{(1)}}+e^{\boldsymbol{\beta}^T\boldsymbol{x}^{(2)}}}$. Note that $1(\delta_1=1)1(y_1<y_2)+1(\delta_2=1)1(y_2<y_1) \leq 1$ and that $w(1-w)\leq0.25$ because $0\leq w \leq 1$. Therefore we have $$\begin{aligned}
    \nu^T\nabla^2_{\beta}\Big\{L^{(s)}(\boldsymbol{\beta})\Big\}\nu&\leq 0.25 \times \nu^T(\bold{x}^{(1)}-\bold{x}^{(2)})(\bold{x}^{(1)}-\bold{x}^{(2)})^T\nu\nonumber\\
    &=0.25\times ||\nu^T(\bold{x}^{(1)}-\bold{x}^{(2)})||_2^2\nonumber\\
    & \leq 0.25\times ||(\bold{x}^{(1)}-\bold{x}^{(2)})||_2^2\times ||\nu||_2^2\nonumber\\
    & \stackrel{(a)}{\leq} 0.25 \times (||\bold{x}^{(1)}||_2^2 + ||\bold{x}^{(2)}||_2^2)\nonumber\\
    & \leq 0.25 \times 2 \times max(||\bold{x}^{(1)}||_2^2, ||\bold{x}^{(2)}||_2^2)\nonumber\\
    & \leq 0.5 \times {\underset{i}{\operatorname{max}}{||\bold{x}^{(i)}||_2^2}},
    \label{B8}\end{aligned}$$ where $(a)$ follows the triangle inequality. Therefore, by assuming our covariates $\boldsymbol{x}$ are bounded, gradient of $L^{(s)}(\boldsymbol{\beta})$ is $D-$Lipschitz-continuous with $D=0.5 \times {\underset{i}{\operatorname{max}}{||\bold{x}^{(i)}||^2}}$. This completes the proof of **Condition 2**. $\square$\
**Proof of condition 3:** The loss function $L(\boldsymbol{\beta})$ belongs to the $C^\infty$ continuous function family and proving is equivalent to proving $$\begin{aligned}
    \exists\, \mu>0, \, s.t.\, \forall \nu \in S_\nu=\{\nu: ||\nu||_2 = 1\},\quad \nu^TI(\boldsymbol{\beta})\nu \geq \mu
    \label{B9}.\end{aligned}$$ where $I(\boldsymbol{\beta}) = \mathbb{E}_{s}[\nabla^2_{\boldsymbol{\beta}}L^{(s)}(\boldsymbol{\beta})]$ is the expected Hessian matrix. Starting from , $I(\boldsymbol{\beta})$ can be written as $$\begin{aligned}
    I(\boldsymbol{\beta}) &=E_s\Big[\nabla^2_{\beta}\Big\{L^{(s)}(\boldsymbol{\beta})\Big\}\Big]\nonumber\\
    &=E_{X, Y, \Delta}\Big[\nabla^2_{\beta}\Big\{L^{(s)}(\boldsymbol{\beta})\Big\}\Big]\nonumber\\
    &\stackrel{(a)}{=}E_{X}\Bigg[E_{Y|X}\Big[E_{\Delta|X, Y}[\nabla^2_{\beta}\Big\{L^{(s)}(\boldsymbol{\beta})\Big\}]\Big]\Bigg]\nonumber\\
    &\stackrel{(b)}{=}(1-p_c)E_{X}\Big[w(1-w)(\bold{x}^{(1)}-\bold{x}^{(2)})^T(\bold{x}^{(1)}-\bold{x}^{(2)})\Big]\nonumber\\
    &=(1-p_c) \int_{X^{(1)}, X^{(2)}}\Big[w(1-w)(\bold{x}^{(1)}-\bold{x}^{(2)})^T(\bold{x}^{(1)}-\bold{x}^{(2)})\Big]P_{X}(x)\nonumber\\
    &\stackrel{(c)}{=} (1-p_c) \int_{\bold{Z}}w(1-w)\bold{Z}^T\bold{Z}P_{\bold{Z}}(\bold{z})
    \label{B10}\end{aligned}$$ where $(a)$ follows the expansion of the intersection using conditional probabilities; $(b)$ follow from the fact that $p_c=E_{\Delta}[1(\delta_i=0))]=P_{\Delta}(\delta_i=0)$ is the probability of censoring and that we have $E_{Y}[1(y_1<y_2))+1(y_2<y_1)] = Pr(y_1<y_2) + Pr(y_2<y_1)=1$; $(c)$ follows from the change of variable $\bold{Z}=\bold{X}^{(1)}-\bold{X}^{(2)}$. Note that since random variables $\bold{X}^{(1)}$ and $\bold{X}^{(2)}$ have density with respect to the Lebesgue measure, random variable $\bold{Z}$ also has a density with respect to the Lebesgue measure. Then we can write $\nu^TI(\boldsymbol{\beta})\nu$ as $$\begin{aligned}
    \nu^TI(\boldsymbol{\beta})\nu &= (1-p_c)\int_{\bold{Z}}w(1-w)\nu^T\bold{Z}^T\bold{Z}\nu P_{\bold{Z}}(\bold{z})\nonumber\\
    &= (1-p_c)\int_{\bold{Z}}w(1-w)||{\nu}^T\bold{z}||_2^2P_{\bold{Z}}(\bold{z})
    \label{B11}\end{aligned}$$

Now we prove strong convexity of $L(\boldsymbol{\beta})$ by contradiction. The negation of statement is: $$\begin{aligned}
    \forall \mu >0,\, \exists \nu_\mu \in S_\nu=\{\nu: ||\nu||_2 = 1\},\, s.t.\quad \nu_\mu^TI(\boldsymbol{\beta})\nu_\mu < \mu.
    \label{B12}\end{aligned}$$

*Claim 1:* Suppose the statement in \[B12\] holds (or equivalently does not hold), then there exists a $\nu^* \in S_\nu=\{\nu: ||\nu||_2 = 1\}$ such that we have $P_\bold{z}(||{\nu^*}^T\bold{z}||_2=0)=1$.\
*Proof:* Since we assumed $\boldsymbol{\beta}^T\boldsymbol{x}$ is bounded, there exists a constant $0<C_w<1$ such that $C_w<w<1-C_w$, and hence $w(1-w)>C_w^2$. Therefore, using , the statement implies that for any $\mu>0$ there exists $\nu_\mu \in S_\nu$ such that $$\begin{aligned}
    \mu &> \nu_\mu^TI(\boldsymbol{\beta})\nu_\mu=(1-p_c)\int_{\bold{Z}}w(1-w)||\nu_\mu^T\bold{z}||_2^2P_{\bold{Z}}(\bold{z})\nonumber\\
    &\stackrel{(a)}{>}(1-p_c)C_w^2\int_{\bold{Z}}||\nu_\mu^T\bold{z}||_2^2P_{\bold{Z}}(\bold{z})\nonumber\\
    &=(1-p_c)C_w^2\int_{\bold{Z}}\Big(1(||\nu_\mu^T\bold{z}||_2^2<\epsilon))+1(||\nu_\mu^T\bold{z}||_2^2>\epsilon)\Big)||\nu_\mu^T\bold{z}||_2^2P_{\bold{Z}}(\bold{z})\nonumber\\
    &\geq (1-p_c)C_w^2\int_{\bold{Z}}1(||\nu_\mu^T\bold{z}||_2^2>\epsilon)||\nu_\mu^T\bold{z}||_2^2P_{\bold{Z}}(\bold{z})\nonumber\\
    &> (1-p_c)C_w^2\epsilon P_{\bold{Z}}(||\nu_\mu^T\bold{z}||_2^2>\epsilon).
    \label{B13}\end{aligned}$$ Note that $\epsilon$ is an arbitrary positive value. Therefore, for such a $\nu_\mu$, expression is equivalent to $$\begin{aligned}
    P_{\bold{Z}}(||\nu_\mu^T\bold{z}||_2^2>\epsilon) < \frac{\mu}{(1-p_c)C_w^2\epsilon}\quad \forall \mu, \, \epsilon>0
    \label{B14}\end{aligned}$$ or equivalently, $$\begin{aligned}
    P_{\bold{Z}}(||\nu_\mu^T\bold{z}||_2^2<\epsilon) >1- \frac{\mu}{(1-p_c)C_w^2\epsilon}\quad \forall \mu, \, \epsilon>0.
    \label{B15}\end{aligned}$$ Thus, there exists an infinite sequence $\nu_1, \dots, \nu_k, \dots$ such that for the choices of $\mu_k(\delta)=\delta(1-p_c)C_w^2\epsilon$ and $\epsilon_k=\frac{1}{k}$ we have $$\begin{aligned}
    \forall\,\delta>0\,, \exists\,K>0\,, s.t.\,, \forall\, k > K: Pr(||\nu_k^T\bold{z}||_2\leq \frac{1}{k}) > 1-\delta.
    \label{B16}\end{aligned}$$ Since $S_\nu$ is a compact space, this sequence has an infinite subsequence converging to a point $\nu^* \in S_\nu$. For such a converging point $\nu^*$, we can write, for all $k$, $$\begin{aligned}
    ||{\nu^*}^T\bold{z}||_2 &= ||(\nu^*-\nu_{k}+\nu_{k})^T\bold{z}||_2\nonumber\\
    & \stackrel{(a)}{\leq} ||(\nu^*-\nu_{k})^T\bold{z}||_2+||\nu_{k}^T\bold{z}||_2\nonumber\\
    & \leq ||\nu^*-\nu_{k}||_2||\bold{z}||_2+||\nu_{k}^T\bold{z}||_2\nonumber\\
    & \stackrel{(b)}{\leq} ||\nu^*-\nu_{k}||_2C_z+||\nu_{k}^T\bold{z}||_2,
    \label{B17}\end{aligned}$$ where $(a)$ follows the triangle inequality and $(b)$ follows the boundedness of variable $\bold{Z}$ (i.e., $||\bold{z}||_2\leq C_z$) due to boundedness of variable $\bold{X}$. From , we may write $$\begin{aligned}
    \forall\,\delta>0\,, \epsilon>0\,, \exists\,K_1>0\,, s.t.\,, \forall\, k > K_1: Pr(||\nu_k^T\bold{z}||_2\leq \frac{\epsilon}{2}) > 1-\delta,
    \label{B18}\end{aligned}$$ and since $\nu_k$ converges to $\nu^*$, we may write $$\begin{aligned}
    \forall\,\delta>0\,, \epsilon>0\,, \exists\,K_2>0\,, s.t.\,, \forall\, k > K_2: Pr(||\nu^*-\nu_k||_2\leq \frac{\epsilon}{2C_z}) > 1-\delta.
    \label{B19}\end{aligned}$$ Then $\forall\,\delta>0$ and $\forall\,\epsilon>0$, for $K_{max}=max(K_1, K_2)$ we have that for all $k>K_{max}$ $$\begin{aligned}
    Pr(||{\nu^*}^T\bold{z}||_2 &\leq ||\nu^*-\nu_{k}||_2C_z+||\nu_{k}^T\bold{z}||_2\nonumber\\
    & \leq \frac{\epsilon}{2C_z}\times C_z + \frac{\epsilon}{2}=\epsilon) > 1-\delta,
    \label{B20}\end{aligned}$$ taking $\epsilon, \delta \rightarrow 0$, we have that $P_{\bold{Z}}(||{\nu^*}^T\bold{z}||_2=0)=1$. This completes the proof of *Claim 1*. $\square$\
*Claim 2:* If random variable $\bold{X}^{(1)}$ and $\bold{X}^{(2)}$ have density with respect to the Lebesgue measure, then we have $P_{\bold{Z}}(||{\nu}^T\bold{z}||_2=0)<1$ for any $\nu \in S_\nu=\{\nu: ||\nu||_2 = 1\}$.\
*Proof:* Random variables $\bold{X}^{(1)}$ and $\bold{X}^{(2)}$ have density with respect to the Lebesgue measure. Therefore, variable $\bold{Z}$ also has a density with respect to the Lebesgue measure. Therefore, for any $\nu \in S_\nu=\{\nu: ||\nu||_2 = 1\}$, $\bold{Z}$ cannot only be on a plane orthogonal to $\nu$. In other words, $$\begin{aligned}
    \forall \nu \in S_\nu=\{\nu: ||\nu||_2 = 1\}, \exists \epsilon_{\nu}>0\, \text{and}\, \exists \delta_{\nu}>0\, s.t.\quad P_{\bold{Z}}(||{\nu^*}^T\bold{z}||_2^2>\epsilon_{\nu}) > \delta_{\nu}.
    \label{B21}\end{aligned}$$ Then for any $\nu \in S_\nu=\{\nu: ||\nu||_2 = 1\}$, the integral $\int_{\bold{Z}}||\nu^T\bold{z}||_2^2dP_{\bold{Z}}(\bold{z})$ can be written as $$\begin{aligned}
    \int_{\bold{Z}}||\nu^T\bold{z}||_2^2dP_{\bold{Z}}(\bold{z})&=\int_{\bold{Z}}\Big(1(||\nu^T\bold{z}||_2^2<\epsilon_{\nu}))+1(||\nu^T\bold{z}||_2^2>\epsilon_{\nu})\Big)||\nu^T\bold{z}||_2^2dP_{\bold{Z}}(\bold{z})\nonumber\\
    &\geq \int_{\bold{Z}}1(||\nu^T\bold{z}||_2^2>\epsilon_{\nu})||\nu^T\bold{z}||_2^2dP_{\bold{Z}}(\bold{z})\nonumber\\ 
    &> \epsilon_{\nu}\int_{\bold{Z}}1(||\nu^T\bold{z}||_2^2>\epsilon_{\nu})dP_{\bold{Z}}(\bold{z})\nonumber\\
    &= \epsilon_{\nu} P_{\bold{Z}}(||{\nu^*}^T\bold{z}||_2^2>\epsilon_{\nu})\nonumber\\
    &> \epsilon_{\nu} \delta_{\nu} > 0,
    \label{B22}\end{aligned}$$ indicating that $P_{\bold{Z}}(||{\nu}^T\bold{z}||_2=0)<1$ (or $P_{\bold{Z}}(||{\nu}^T\bold{z}||_2=0)\not=1$). This completes the proof of *Claim 2*. $\square$\
*Claim 1* indicated that if the strong convexity statement does not hold, we must have $P_\bold{z}(||{\nu^*}^T\bold{z}||_2=0)=1$. This is in contradiction with the result of *Claim 2* indicating that for random variable $\bold{Z}$ with Lebesgue density we have $P_\bold{z}(||{\nu^*}^T\bold{z}||_2=0)\not=1$. Therefore, is not true and the statement holds, i.e., $L(\boldsymbol{\beta})$ is strongly convex. This completes the proof of **Condition 3**. $\square$\
**Proof of condition 4:** The gradient of $L^{(s)}(\boldsymbol{\beta})$ may be written as $$\begin{aligned}
    \nabla_{\boldsymbol{\beta}}L^{(s)}(\boldsymbol{\beta}) = \sum_{i=1}^s1(\delta_i=1)(\bold{x}^{(i)}-\sum_{j \in R_i}w_j\boldsymbol{x}^{(j)}).
    \label{B23}\end{aligned}$$ Then we can write $$\begin{aligned}
    ||\nabla_{\boldsymbol{\beta}}L^{(s)}(\boldsymbol{\beta})||^2 &= ||\sum_{i=1}^s1(\delta_i=1)(\bold{x}^{(i)}-\sum_{j \in R_i}w_j\boldsymbol{x}^{(j)})||^2 \nonumber\\
    &\stackrel{(a)}{\leq} \Big(\sum_{i=1}^s(||1(\delta_i=1)||\times||\bold{x}^{(i)}-\sum_{j \in R_i}w_j\boldsymbol{x}^{(j)})||\Big)^2\nonumber\\
    &\stackrel{(b)}{\leq} \Big(\sum_{i=1}^s(||\bold{x}^{(i)}||+||\sum_{j \in R_i}w_j\boldsymbol{x}^{(j)}||)\Big)^2\nonumber\\
    &\stackrel{(c)}{\leq} \Big(\sum_{i=1}^s(||\bold{x}^{(i)}||+ {\underset{j \in R_i}{\operatorname{max}}{||\boldsymbol{x}^{(j)}||)}} \Big)^2\nonumber\\
    &\leq 4S^2{\underset{i}{\operatorname{max}}{||\boldsymbol{x}^{(i)}||^2}}
    \label{B24}\end{aligned}$$ where $(a)$ follows from the triangle inequality; $(b)$ and $(c)$ follow from the triangle inequality and that $w_j\leq1$. Therefore, given boundedness of covariates $\bold{x}$, we can choose $\sigma=2S{\underset{i}{\operatorname{max}}{||\boldsymbol{x}^{(i)}}}||$. This completes the proof of **Condition 4**. $\square$\
We showed that all four conditions are satisfied and thus the results in [@Eric2011] give us our claimed convergence rates for both single SGD-based estimates and averaging over estimates (Polyak-Ruppert average).

Implementation of streaming and non-streaming algorithms using BigSurvSGD {#appC}
=========================================================================

In this section, we present the implementation of both streaming and non-streaming mini-batch stochastic gradient descent algorithms using our proposed framework BigSurvSGD. Without loss of generality, we only present algorithms without the moment-based step-size adaptation.

Implementation of a streaming algorithm {#secAC1}
---------------------------------------

As we discussed before, our proposed framework facilitates the implementation of an algorithm dealing with the streaming data. Such an implementation is very straightforward. The coefficient estimate (i.e., $\boldsymbol{\beta}$) is updated in a streaming fashion. Thus there is no need to collect all the data before estimation (as would be required by [coxph()]{}. Each mini-batch can be as small as the size of the strata we use. Thus we will not run into memory issues. Algorithm \[algStreaming\] gives the details of a streaming algorithm using our proposed framework.

**Initialization**:\
Choose strata size $s$\
$\hat{\boldsymbol{\beta}}(0)=\boldsymbol{0}$\
Choose $C$ for $\alpha_m=\frac{C}{\sqrt{m}}$

Implementation of a non-streaming mini-batch gradient descent algorithm {#secAC2}
-----------------------------------------------------------------------

In many cases, it is of interest to consider each datum more than once (this can empirically improve performance). In this case, we will still use a mini-batch stochastic gradient descent algorithm. Here we present the implementation of a non-streaming mini-batch gradient descent algorithm (that includes mini-batches and multiple epochs). We being by splitting our data evenly into mini-batches, each with $K$ strata of size $s$. Then we iteratively update the estimate using batches of strata. Using multiple, rather than single strata per batch results in more stable (less noisy) gradients. We use the learning rate $\frac{\alpha_m}{s\times K}$ instead of $\frac{\alpha_m}{s}$ as the gradient in each step of mini-batch gradient descent is the sum of gradients over $K$ strata with size $s$. Algorithm \[algMiniBatch\] summarizes the implementation of a non-streaming mini-batch gradient descent algorithm using our proposed framework. One strength of such an implementation compared to [coxph()]{} is we can read batches of strata chunk-by-chunk from the hard drive. Therefore, such an implementation can handle large amounts of data without facing memory issues.

**Initialization**:\
Choose strata size $s$\
Choose batch size $K$\
Choose number of epochs $n_E$\
$m=0$\
$\hat{\boldsymbol{\beta}}(0)=\boldsymbol{0}$\
Choose $C$ for $\alpha_m=\frac{C}{\sqrt{m}}$\

0.2in
