FastSurvival: Hidden Computational Blessings in Training Cox Proportional Hazards Models
Authors: Jiachang Liu
Rui Zhang
Cynthia Rudin
Venue: NeurIPS 2024 poster (2024)
Keywords: survival analysis
cox proportional hazards model
optimization
first-order method
sparse learning
Primary area: interpretability_and_explainability

Abstract: Survival analysis is an important research topic with applications in healthcare, business, and manufacturing. One essential tool in this area is the Cox proportional hazards (CPH) model, which is widely used for its interpretability, flexibility, and predictive performance. However, for modern data science challenges such as high dimensionality (both $n$ and $p$) and high feature correlations, current algorithms to train the CPH model have drawbacks, preventing us from using the CPH model at its full potential. The root cause is that the current algorithms, based on the Newton method, have trouble converging due to vanishing second order derivatives when outside the local region of the minimizer. To circumvent this problem, we propose new optimization methods by constructing and minimizing surrogate functions that exploit hidden mathematical structures of the CPH model. Our new methods are easy to implement and ensure monotonic loss decrease and global convergence. Empirically, we verify the computational efficiency of our methods. As a direct application, we show how our optimization methods can be used to solve the cardinality-constrained CPH problem, producing very sparse high-quality models that were not previously practical to construct. We list several extensions that our breakthrough enables, including optimization opportunities, theoretical questions on CPH's mathematical structure, as well as other CPH-related applications.

FastSurvival: Hidden Computational Blessings in Training Cox Proportional Hazards Models
Jiachang Liu', Rui Zhang", Cynthia Rudin
CComell University, Duke University
j jiachang www.lllinnee..... edu, x z.aang@dukee.co edu, cynthia@cs, duke. edu
Abstract
Survival analysis an important research topic with applications in healthcare, business, and manufacturing. One essential tool in this area is the Cox proportional hazards (CPH) model, which is widely used for its interpretability, flexibility, and predictive performance. However, for modern data science challenges such as high dimensionality (both and p) and high feature correlations, current algorithms to train the CPH model have drawbacks, preventing us from using the CPH model its full potential. The root cause is that the current algorithms, based on the Newton method, have trouble converging due to vanishing second order derivatives when outside the local region of the minimizer. To circumvent this problem, we propose new optimization methods by constructing and minimizing surrogate functions that exploit hidden mathematical structures of the CPH model. Our new methods are easy to implement and ensure monotonic loss decrease and global convergence. Empirically, we verify the computational efficiency of our methods. As direct application, we show how our optimization methods can be used to solve the cardinality-conster constraiee CPH problem, producing very sparse high-quality models that were not previously practical to construct. We list several extensions that our breakthrough enables, including optimization opportunities, theoretical questions on CPH's mathematical structure, as well as other CPH-related applications.
Introduction
Survival analysis, which studies time-to-event data, is an important research topic with a wide range of real-world applications. In medicine, survival analysis has been employed to model when patient will die 444 38 [0]. business, it useful for attrition prediction (30) (when an employees resigns) and churn prediction (31 (when a customer unsubscribes), and in manufacturing, it is used to predict when a physical system breaks down 42 50]. fundamental tool analyzing such data the Cox
proportional hazards (CPH) model 8 linear model under the assumption that features have multiplicative effect on the risk of failure/event. Simple yet powerful, the CPH model has enjoyed great popularity due to its modeling flexibility (when coupled with additive models 27 28 Moreover, in contrast to black box models, is both interpretable and accurate.
However, with the advent of larger sample and feature spaces and more complex data, new challenges arise using the CPH model to its full potential. Ideally, practitioners want to produce CPH models repeatedly, with feature engineering and preprocessing between iterations. Additionally, they want CPH model to identify important variables (64 [4], even in presence of highly correlated features. However, current optimization methods for training the CPH model do not meet these needs. Current
algorithms 662 22 23 54), based on the generic Newton's 's method, are computationally intensive.
More importantly, due to both vanishing second-order derivatives [53] and the use of approximation strategies that trade precision for efficiency, existing optimization methods have trouble converging. either with the loss blowing up or the algorithm converging very slowly when require the precision
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

necessary to handle correlated variables. The latter issue is the core reason for incorrect variable selection when features are highly correlated.
In this work, we propose new optimization methods to train the CPH model and show that there is not necessarily prefision-effcieency tradeoff. Despite the CPH model being seemingly amenable to classical optimization approaches such as coordinate descent, this has not been attempted for the original CPH loss function due its daunting complexity. However, through careful examination, we show instead that the complexity of the CPH loss function is really a blessing, rather than a curse. We discover hidden mathematical structures that allow us to design very efficient algorithms. We show both the first and second- order derivatives at each coordinate can be computed exactly in linear time complexity (7)). Moreover, we show both derivatives are Lipschitz- cotthitz--ontuuuuss by making novel connections with the second and third central moment calculation in probability theory and statistics. All these discoveries lead us to design algorithms that essentially minimize quadratie surrogate function and a cubic surrogate function, respectively. They are extremely easy to implement.
Empirically, we demonstrate the superior speed of our algorithms on large-scale datasets. In general, ours are significantly faster than all existing methods and rapidly converge to optimal high-precision solutions. Because our methods produce high-quality solutions, we apply them for variable selection in challenging regimes where features are highly correlated. We solve difficult cardinalit--constrannee constraine CPH problems and produce models that are much sparser than the state-of-the ftt--hhettttt art methods.
In summary, our contributions are: (1) We find a critical flaw in the current optimization algorithms for the CPH method by pinpointing that they converge slowly with low precision. Sometimes, the loss does not converge and explodes. (2) To circumvent this issue, we propose novel algorithms that minimize quadratic and cubic surrogate function, respectively, with guaranteed convergence and loss descent at each iteration. The core novelty lies in discovering hidden mathematical structure, which allows for an efficiency way (O))) of calculating the second-order partial derivatives exactly. In addition, we show the first and second order partial derivatives are Lipschitz continuous. To calculate these Lipschitz constants, we leverage second and third central moments from theoretical statistics and probability theory. (3) Empirically, our method enjoys fast speed in training the loss function and results in superior performance when solving ceedinallty-consrramnn problems.
Our work constitutes a methodological breakthrough in training CPH models. At the end of the paper, we also discuss several exciting extensions and follow-up questions, showing how our new perspectives and discoveries open doors to many new research opportunities.
Preliminaries
Given a time-to-event dataset of n samples with (x,, t,, 8 where x E R' is the feature vector with length p, TR. is the observation time, and 5, {0, 1] is indicator with 1 indicating that a failure event has happened, the CPH model can be used to learn and predict the risk of failure, commonly known the hazard function in survival analysis. The CPH model predicts the hazard (t()) for sample i in a semiparametric way 622. For review of related work, please see Appendix B
ht(1) ho (0)e
where An (4) is a baseline hazard function shared by all samples, and 3 E IRP is the parameter of interest. The nice thing about the CPH model is that (4) cancels out if we look only at the ratio of hazards of sample i vs. all remaining samples at time t, ie..
hi(ts) Ejen, (y(44)
c
ΣSER, "
where R; = [j 1j 2 43 is the set of indices whose observation time is greater than equal that of sample i. Such ratios are also called partial likelihoods. To estimate the parameter of interest, P maximize the joint partial likelihood of all samples with failure events, which can be written as
T
L(Ø) = 1116,
(3)
Σ EjER, "
2

---

This is equivalent to minimizing the negative log partial likelihood 62] which is defined as
((B) = - log L(3) = Σ d2 log
Σ
c 8
T 3
(4)
jER,
The loss function ((3), while convex, is very mathematically involved. In addition to the double sum, the inner summation over j is with respect to a different index set R2 for each outer summation index i. Such daunting complexity makes it difficult to employ first-order methods such as gradient descent because we canmot easily pick the right step size for each iteration, which plays a crucial role in practical running time. Therefore, past efforts have been focused on developing Newton-type (second-order) methods, where the loss function is approximated by second-order Taylor expansion:
[(3 + A3) ~ {(3) + V)((3)) +
1
== f(43).
(5)
The function f(4,3) can be minimized by solving a linear system: A33 = - To reveal the computational nuances more explicitly, we use an intermediate variable 77 with 77 X/3. Then we can rewrite the approximation function (A3 as:
((AB) {B + A3) ((n) + v,,(())² XA3 + 2 AB (m)
(6)
At each iteration, caleulating the Hessian matrix (n))(n) requires O(n²) complexity. Past methods on the CPH model have resorted to various approximation strategics by replacing (2,(()) with (()) to reduce the computational complexity:
1. Exact Newton
(nn) = V (n)
# no approximation
[ (n))a
if otherwise
2. Quasi Newton
H(n);; 0
# ignore off diagonal terms
3. Proximal Newton
H()) = diag(V vv, 8),
# diagonal upper bound on V (())
where diag(-) constructs matrix with its diagonal equal to the imput vector and other entries equal to There are two major problems with the above approaches. One common problem is that the these Newton-type methods inherently have trouble converging beyond the local region of minimizers without backtrack line search [53]. We provide a concrete example to demonstrate this issue in the experiment section. Ideally, we want to avoid backtracking because this increases the running time. In contrast, our methods do not have this flaw and guarantee global convergence.
The other problem is that when the above approaches do converge to the optimal solutions, none of them can converge with high precision fast enough (in a practical sense). The exact Newton's s method 222 has local quadratic convergence rate, but each iteration can take a long time. Quasi Newton 62 and proximal Newton [51] methods are computationally much cheaper to evaluate per iteration, but they make less progress toward the optimal solution. In the next section, we show that, by exploiting hidden mathematical structure, we can obtain the best of both worlds: cheap evaluation per iteration and fast convergence with respect to the number of iterations.
3 Methodology
3.1 New Formulas for First, Second, and Third Order Partial Derivatives
As we have mentioned, the reason for the diagonal approximations of (Z(((n) is to reduce the complexity of the mathematics and associated high computational cost. Here, we take a completely different approach from past methods. First, we avoid making any approximations and embrace the full Hessian matrix. Second, we bypass the intermediate step of calculating the Hessian in the sample space 7) and focus on the Hessian in the feature space B. The involved mathematics may already sound complicated, but we do not stop here. We apply these two ideas not only to the second order partial derivatives but the third order partial derivatives as well. Although this seems like burdensome task, we show that the end result very elegant and has an intuitive interpretation. Out of complexity comes simplicity. We summarize the relevant results in the first following theorem:
'See the skglm tutorial at https://contrrb.sccittt.t learn, http:///ww...//.kk///////o/////a ittaff
Htaths cox-dataafi

---

Theorem 3.1. For the CPH loss function defined in Equation 4). partial derivatives with respect coordinate are: 1st order partial derivative:
the first, second, and third order
0E(33
e
Σd
Σ Rr
Xxi
Σ
6,
08t
ZJER,
e
i=1
WE
i=1
2nd order partial derivative:
)2((3)
Th
e X
Σd
X
Σ
(8)
a
EjEn, e
Ejen, eni
KER,
KER
3rd order partial derivative:
3
6
Σ
eh
x3 Xki+2 +2
Σ
c
XAil
081
EJER, e
E
jjER e
==1
LEE;
MER
-3 Ejen, ew
e
X
Σ NER
Enen
e
Xm
(9)
Ejen, e
KER,
The proof can be found in Appendix A The first, second, and third order partial derivatives all have probabilistic interpretation. Notice that for any i, the coefficients in front of XKt. X2, and X2, nonnegative and sum up to 1, i.e., e /(>ERR, e":) 2 and ΣkER, R. [e /(jeR, e*))] = 1. Then, we can regard these coefficients as a discrete probability distribution. Thus, for Equation
the term inside [-] resembles the variance or second order central moment formula: E[X (B[X) )
ET(X - [] For Equation (9), the term inside resembles the skewness or third order central
moment formula: [[X*] + 2[E[X)) 3 [[[( E[X])*]
One may wonder whether for higher orders (order r 2 4), the relationship between the r-th order partial derivative and r-th central moment still preserve. The answer is no and this can be easily deduced from the following lemma. The proof can be found in Appendix
Lemma 3.2. Let us define Cr to be the r-th central moment with
C= Σ
th
Σ
t
XmL1
(10)
EjER, e
EjER, e
RER.
RIER.
Then we can calculate the partial derivative of Cr with respect to B1 as:
i 03, Cr Cr+1 - - -0- Cr-1-
(11)
From Lemma 3.2 we can see why the connection to central moment does not work for higher order partial derivatives. Therefore, we get /0, C3. However, for 3, C,-1 in general is not zero, so we cannot extrapolate this pattern to higher order partial derivatives.
2, the second term in Equation (TT) disappears, i.e., C,-1 C1
Theorem 3.1 forms the basis upon which we build everything else. These results are not only mathematically interesting but also have significant implications for computation, which we elaborate in the next two sections.
3.2 Time Complexity of First and Second Order Partial Derivative Calculation
From the connections to the second and third central moment, have the following corollary regarding the time complexity of calculating the first, second, and third order derivatives:
Corollary 3.3. For the CPH model, the time complexities to calculate
0((3) 00:
and ('
are O(n).
age

---

This is surprising result, especially for the second order partial derivatives. The intermediate Hessian, V (n)) takes (n²²) to compute, so we would expect the second order partial derivative, (2 = ej x (()))(()) would take O(n²) to compute as well. Yet. the time complexity is 88; just O(n). We use the first order partial derivative formula, Equation 77 as an example to explain why this happens. We ignore the second term Di= 6, because it is just constant. The first term in Equation 7 can be rewritten as:
> di EMER e Xkrl 1=1
Σ RIER,
LJER, eb Xky!
(12)
ΣJER, C²
I=1
Let us focus on the numerator inside the parenthesis for now. For the entire sequence (i = 2..,.... m) of numerator terms, cumulative summation. The same is true when we obtain the entire sequence of denominators. Once we have all numerators and denominators. calculating the entire sequence of ratios inside the parenthesis also costs Onn). Finally, multiplying each ratio with ó; and summing up all these products costs O(n) well. Therefore, the computational cost calculate the first order partial derivative is O(n). We can apply the same idea to the second order partial derivative formula, Equation & to show that the computational complexity is also (n)). Note that the reverse cumulative summation trick has already been explored in [62] for calculating the diagonal of V (() in the sample space v), but this trick has not been used calculate the partial derivatives in the feature space B.
obtain all of them together at the cost of ((n) by performing reverse
We will later see how this (nn) time complexity allows us to design a second order optimization method whose evaluation cost per iteration is just cheap first order optimization method. Before we discuss that, let us continue and discuss another computational implication of Theorem 3.1
3.3 Hipschitz-Continuigy Property of First and Second Order Partial Derivatives
The connection to the central moment calculation allows us to conclude that the first and second order partial derivatives are Lipschitz-continuuus
Moreover, we can calculate these Lipschitz constants explicitly. Recall that for a univariate function fz)), we say that the function is Lipschitz- continuous (5] there exists Z such that for any two points in the domain of f, i.e., x, y E D()), we have [(x)) [w) < L|3: y/- The value L is called the Lipschitz constant for function f()). If the function is continuously differentiable, the previous definition is equivalent to the condition where the first order derivative is bounded, [x) s L for any E R (5).
Not only can we say that the first and second order partial derivatives are Lipschitz-continuov but we can also calculate the Lipschitz constants explicitly. We summarize the results in the theorem below. The proof can be found in Appendix A
Theorem 3.4. For the second order partial derivatives in Equation bounded by the following formula:
its absolute values are
0
Σ
5,
max RIER.
Xx1
min k.ER.
Xx,1)²
(13)
037
i=1
For the third order partial derivatives in Equation formula:
its absolute values are bounded by the following
(Pe(m)
1 max Xb, KjER;
3
Σ
8
min Xk,1
(14)
033
6v
R
The availability of these Lipschitz constants suggests that we might construct surrogate functions.
3.4 Quadratic and Cubic Surrogate Functions
We now have all the tools at hand to attack the original optimization problem. For a univariate convex fx), if we have access to L2, the Lipschitz constant for its first order derivative, then can construct the following quadratic surrogate function (-) [53]:
1
f(x + Ar) S f())) f Ar + L2Ax² =: =.
(15)
2
5

---

If we have access to L3. the Lipschitz constant for its second order derivative, then we can construct the following cubic surrogate function h. (-) [33]
f(x +x) S ft + ")" " -
h (Ax)
(16)
nice thing about these surrogate functions is that their minimizers can be computed analytically:
argmin Iu (Ax) - (((x)
(17)
Az
argmin hy (Ax) = sgn(f(x))) (= ((f + 2L
(18)
L3
where the function sgn(.) extracts the sign (+ or -) of the input. The analytical solution to the quadratic surrogate function is well known, but the analytical solution this cubic surrogate function has not been well studied. We provide a derivation for Equation (18) ( in Appendix [A]
Since these surrogate functions are convex and are upper bounds of the original functions, minimizing them will lead to a decrease of the original function fxx as well. This explains why our methods ensure monotonic decrease in loss and guarantee global convergence. The final algorithms are very easy to understand and can be thought as coordinate descent etppentttypeeaeee methods. We anticipate these core ideas be applied solve wide range of problems related the CPH model. In the next subsection, we showcase two problems our algorithms can tackle.
3.5 Applications to Regularized and Constrained Problems
Regularized Problem The first problem is the regularized CPH problem whose penalty terms are separable. The penalties that qualify for this category include LASSO [64]. ElasticNet [72] SCAD [5I, MCP (68). etc. For the e-regularized problems, we can fact find analytical solutions
For the quadratic surrogate function, solving the fi-regularized problem Equation (15) is equivalent solving the following optimization problem (with = f'(x), Lz, and c 2:),
1 bAr + +1lc+ Arl.
A argmin aAx +
(19)
Az
2
The solution for the above problem is
\1//b if be a < - (a
A#
(a A1)/b bc- X1
(20)
otherwise.
For the cubic surrogate function, solving the f1-regularized problem of Equation (66 is equivalent to solving the following optimization problem (with a = f'()), f"())) ( L3. and = x):
1 bb2² 6
A argmin aAz +
dl/1 A Ar/,
(21)
Ax
2
whose solution is:
sgn (d)
- b + 1 (sgn( (d)a + \1)
/c if sgn(d)a + \1 5 0
sgn ())
6+ / 62 22( d)a A1)
/c sgn(d)(a bd) Jo > A
A
(22)
sgn d)
b / 82 + 20(sgn d)a X1 /c if sgn(d)(a bd) ed -
otherwise.
Equation 20 is well known a slightly different format. Equation 22 has not been well studied in the past. We provide derivations for both in Appendix A
For the ElasticNet problem where the penalty is A A1 + 2211-1 we can also easily obtain analytical solutions. The trick is to absorb the first and second order derivatives into the coefficients of the surrogate functions and only have Au|I-I| penalty.

---

Constrained Problem The second problem is the cetiiiallty-oonstaiinnn CPH problem. Recently, the beam search framework combination of the beam search method [66] from natural language processing and generalized orthogonal matching pursuit [13) has shown promise in finding near opti- mal solutions for a class of lg-constrained nonconvex problems, including sparse ridge regression [47
and sparse logistic regression [48]
Similar to the generalized orthogonal matching pursuit algorithm, we expand our support (starting from an empty set) by adding one feature at a time until the cardinality is satisfied. However, instead of selecting features based on partial derivatives, we select features based on which coefficient, if optimized, can result in the largest decrease of the loss function. After the feature is added into the support, we fine-tune all nonzero coefficients in the support. Additionally, during each support expansion step, we select multiple feature candidate instead of the best one, similar to the core idea in beam search. We use our coordinate descent methods to solve the feature selection step and the coefficient fine-tuning step.
Although the beam search framework has already been proposed for other aardinality-constraine problems, it cannot be applied directly to the CPH model without our coordinate descent methods to select features, especially in the highly correlated settings.
4 Experiments
We test the effectiveness of our optimization methods on both synthetic and real-world datasets. We experiments for both regularized and constrained problems mentioned in Section (.5 Our main objectives are: 1) When minimizing the same objective functions, how fast can our methods converge to the optimal solutions when compared with all existing optimization methods for the CPH model? 2) When coupled with the beam search framework, how well can our methods help with variable selection when compared with the state-of of-the-arts methods, especially for challenging scenarios where features are highly correlated?
4.1 Accessing How Fast Our Methods Converge to Optimal Solutions
We compare our methods (one based on the quadratic surrogate function and the other based on the cubic surrogate function) with the existing optimization methods outlined in Section [ exact Newton method, the quasi Newton method, and the proximal Newton method. We run on both fz-regularized CPH problems and & + Ez-regularized CPH problems. The choices of these regularizations are: A = (0, }} and X₂ = 51. The coefficients are all initialized to be O. In the main paper, we show results on the Flchain dataset in Figure ! More results on other datasets can be found in Appendix D) During each iteration, the bascline methods 162 51 optimize all coefficients at once, whereas methods optimize coefficients sequentially with respect to the original loss function. To assess the per-iteration convergence rate, we plot the CPH loss against the number of iterations. To assess the practical running speed, we plot the CPH loss against the overall time elapsed (wall clock). From the left two plots of loss vs. number of iterations, we see that the Newton-rype baselines sometimes have losses that blow up or increase during the initial phase of optimization. This is a common problem of Newton's 's method. Our methods are the only ones with monotonically decreasing loss curves. This is the main reason why only our method can be used for the beam search framework in the variable selection experiments. From the right two plots of loss overall time elapsed, methods are significantly faster than the baselines. This is due to the fact that both our first and second order partial derivatives are very cheap to compute (with time complexity (nz)), as we have explained in Section 3.2
that our
4.2 Accessing How Well Our Methods Perform Variable Selection
We compare our method with both Cox-based methods and other model classes. For Cox- baseb models, we run on both synthetic datasets and real-world datasets. For other model classes, we only run on the real-world datasets. To assess how well different methods select important variables, features are highly correlated all datasets. Synthetic datasets are generated with high correlation level, p = 0.9. For each continuous feature on the real-world datasets, we perform binary thresholding for preprocessing [49] to obtain many one-hot encoded binary features. This preprocessing step result in highly correlated features on which it is challenging to perform variable selection. We use the
7

---

Loss Iters
Loss Time
Loss Iters
Loss Time
lambdda lambda2=1
lambdel lambda2=1
lambdal lambda2=5
lambdal 11, lambda2-5
15003
15000
16000
14000
losses
15500
11500
13000
13000
12300
12500
12000
12000
10
10f
10
0.0
255
5.0
2.5
10.0 1255
iter's
elapsed times (sec)
Iters
ela peopsed_times (sec)
NewtonExect
NewtonProx
NewtonQuasi
QuedraticSunogate (Ours)
CubicSurrogtee (Ours)
Figure 1: Efficiency experiments on the first fold of the Flchain dataset. a) The left two plots are on the fz-regularized problem with 22 = For all Newton-type methods, the losses blow up when regularization is weak. In contrast, our methods (quadratic and cubic surrogates) ensure monotonic decrease of losses. b) The right two plots are on the & + fz-regularized problem with ' = 1 and A2 5. The exact Newton method cannot be directly applied, so we compare only with quasi Newton [62] and proximal Newton [51] methods, which have losses that increase at the beginning. Our methods are significantly faster than both baselines. Because the evaluation cost per iteration is very cheap for our methods, we are significantly faster in terms of wall clock time (see the difference between the third and fourth plots). See Appendix D for results other datasets.
Support Size vs fl score Dataset: SyntheticHinnCrrHihhhmmm
Support Size vs fl score
Support Size vs fl score Datasel: SyntheticHighCorHighhho
Datasel: SymtheticHgghorrHighhh
1.0
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
10a
10a
10a
10a Support Size
10a
10a Support Size
Support Size
skgimAlassoooo
abess
sesurvCaxnet
ReamSearch (Duns)
Figure 2: Variable selection on synthetic datasets with high correlation (correlation level p p= 0.9). From left right, the sample sizes are 1200, 1000, and 800, respectively. The Fl score (the higher the better) is closely related to the support recovery rate, On the left two plots, we can see our method recovers the true variables significantly better than other methods (100% recovery rate on the left plot; true support size is 15). As the sample size decreases, the Fl score decreases for all methods.
following metrics to evaluate our solution qualities: CPH loss, CIndex, and IBS. On the synthetic datasets where we know the true coefficients, we also calculate the Fl score. We perform 5-fold cross validation and report the mean and standard deviation of different metrics on both the training and test sets. For details about the experimental setup, please see Appendix C For Cox-based methods, we compare our method with Coxnet, Abess, and Adaptive Lasso.
Results on the synthetic datasets are shown in Figure 2 We plot support size vs F1 score. The F1 score is closely related to the support recovery rate. Our method performs significantly better than the basclines. In particular, on the leftmost plot with 1200 samples, our method is the only one to achieve 100% recovery rate; the true support size is 15 and recover all 15 features with a model of size 15. Results for the Employee Attrition dataset are shown in Figure 3 We plot support size vs. Cindex and support size vs. IBS. Similar to the trend on the synthetic datasets, our method performs significantly better than the baselines in terms of both metrics. Lastly, we compare our method with other model classes on the Dialysis dataset. The results are shown in Figure 4 We plot support size vs. CIndex and support size vs. IBS. The results indicate that other model classes are prone to overfitting on the training sets. Our method achieves the best eecuracy-sparsic tradeoff. We are able to obtain solutions with the smallest number of coefficients without losing predictive performance.
8

---

Support Size Vs trein, cindex Dataset: Emplyyssrrtiion
Support Size vs test cindex Dataset: Employerátrition
Support Size vs train ibs Dataset: Employecdtrrtian
Support Size vs test lbs Dataset: Employecstriition
0.90
0.00
0.18
0.80
0.88
ttain cindex
0.60
cindex
0.86
0.36
0.84
0.14
10.02
0.13
0.00
6:10
0.10
0.78
0.76
20 Support Size
10
30
20
10
10
0.00
10
70
30
Support Size
Support Size
Support Size
seglin ALassooo
abiss
SksurCCesret
BeamSearchh (Durs)
Figure 3: Variable selection on the Employec Attrition dataset. We show support size vs. CIndex (left two plots, the higher the better) and support size vs. IBS score (right two plots, the lower the better). We compare our method with Cox- based sparse learning methods. For both metrics, our method is significantly better than other baselines.
Support Size train cindex
Support test cindex Dataset: Dialysis
Support Size Dataset: Dialysis
Support Size test ibs Dataset: Dialysis
Dataset: Dialysis
1.10
0.95
0.20
Vain clidex
0.50
0.19
0.85
0.125
0.60
0.00
0.100
0.17
0.75
0.075
0.60
200000 600000 600000
0.8
103
101
10 Support Size
Support Size
Support Size
Support Size
iksurviree
sksurvii
stsurvhinvee00
sksurviasttiV
BearSSearch (Guns)
Figure 4: Variable selection on the Dialysis dataset. We show support size vs. Cindex (left two plots, the higher the better) and support size vs. IBS score (right two plots, the lower the better). We compare our method with other model classes. For both metrics, our method obtains solutions that are significantly sparser than other model classes without losing accuracy on the test sets. Other model classes are prone to overfitting on the training sets.
All these results demonstrate the superior sparse learning capability of our method. For more results, with all baselines on all datasets, please see Appendix D
Limitations of FastSurvival Our work focuses on efficient training and effective variable selection of the CPH model. Other model classes, such as trees, random forests, and neural networks, have their own unique merits in capturing complex patterns when the linear (or in our case, additive) model assumption is not satisfied. Another limitation using the CPH model itself, since its assumptions do not always hold. Handling this question is out of scope for this work.
Conclusion and Future Outlook
We presented new optimization methods train the Cox proportional hazards (CPH) model by con- structing and minimizing either a quadratic or a cubic surrogate function. We achieve computational efficiency by exploiting the hidden mathematical structures discovered for the CPH model. Our algorithms are able to train the model significantly faster than previous approaches while avoiding the issue of loss explosion. Furthermore, when applied to the variable selection problem, our method can produce solutions with much fewer parameters while maintaining predictive performance. There are many possible extensions to build upon this work. On the optimization side, will be interest- ing to see whether we can derive analytical solutions for other types of regularizers mentioned in Section 3.5 On the theoretical side, questions remain whether higher order partial derivatives are Lipschitz-continuuus and how to compute these Lipschitz constants. On the application side, we can apply our method to solve the CPH models with time-varying features [6 stratifications (40, and feature interactions [45].

---

Code Availability
Implementations of FastSurvival discussed in this paper are available https:// com/ jjachangjln//aattturrvvaa
Acknowledgements
The authors gratefully acknowledge funding support from the National Institutes of Health under 5R01-DA054994 and the Department of Energy under DE-SC0021358.
References
[1] A. Bender, A. Groll, and Scheipl. generalized additive model approach to time-to-event analysis. Statistical Modelling, 1893-4):29--32,, 2018.
[2] Q. Bertrand, Q. Klopfenstein, P.-A. Bannier, G. Gidel, and M. Massias. Beyond 11: Faster and better sparse models with skglm. Advances in Neural Information Processing Systems, 35:38950-38965, 2022.
[3] D. Bertsimas, J. Dunn, E. Gibson, and A. Orfanoudaki. Optimal survival trees. Machine Learning, 11 (11())2951-3022 2022.
[4] D. Bertsimas, A. King, and R. Mazumder. Best subset selection via a modern optimization lens. The Annals of Statistics, pages 813-852, 2016.
[5] S. Bubeck al. Convex optimization: Algorithms and complexity. Foundations and Trends@ in Machine Learning, 8(3-4):231-357 2015.
[6] Z. Che, S. Purushotham, K. Cho, D. Sontag, and Y. Liu. Recurrent neural networks for multivariate time series with missing values. Scientific Reports, 8(1):6085, 2018.
[7] T. Ching, X. Zhu, and X. Garmire. Cox-nnet artificial neural network method for prognosis prediction of high-throughput omics data. PLoS Computational Biology. 14(4):e100607) 1006076, 2018.
[8] D. R. Cox. Regression models and tables. Journal of the Royal Statistical Sociery: Series (Methodological), 34(2):187-202, 1972.
[9] I. Daubechies, M. Defrise, and C. De Mol. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 57(1 (((1)): 11:1431-1457 2004.
[10] A. Dispenzieri, J. A. Katzmann, R. A. Kyle, D. R. Larson, T. M. Themeau, C. L. Colby, R. J. Clark, G. P. Mead, S. Kumar, J. Melton III, et al. Use of nonclonal serum immunoglobllin free light chains predict overall survival in the general population. In Mayo Clinic Proceedings, volume 87, pages 517-523. Elsevier, 2012.
[11] V. B. Djeundje and J. Crook. Identifying hidden patterns in credit risk survival data using generalised additive models. European Journal of Operational Research, 277(1):366-36,, 2019.
[12] E. Drysdale. Survset: An open-source time-to-event dataset repository. arXiv preprint arXiv:2203.0309. 2022.
[13] E. R. Elenberg, Khanna, A. G. Dimakis, and S. Negahban. Restricted strong convexity implies weak submodularity. The Annals of Statistics, 46(66):3539-33)) 3568, 2018.
[14] J. Fan, Y. Feng, and Y. Wu. High- nigh-dimenssonal variable selection for Cox's 's proportional hazards model. In Borrowing Strength: Theory Powering Applications-A Festschrifi for Lawrence Brown, volume 6, pages 70-87. Institute of Mathematical Statistics, 2010.
[15] J. Fan and R. Li. Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American statistical Association, 1348-1336)) 2001.
[16] I.. D. Fisher and Y. Lin. Time-dependent covariates in the cox proportional-hazaree regression model. Annual review of public health, 20(1):145-1577 1999.
[17] S. Fotso et al. PySurvival: Open source package for survival analysis modeling. 2019.
10

---

[18] J. Friedman, T. Hastie, H. Höfling, and R. Tibshirani. Pathwise coordinate optimization. The Annals of Applied Statistics, 1(2):302 332, 2007.
[19] Friedman, Hastie, and R. Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 933 2008.
[20] J. Friedman, T. Hastie, and R. Tibshirani. Regularization paths for generalized lincar models via coordinate descent. Journal of statistical software, 33(1):1, 2010.
[21] E. Giunchiglia, A. Nemchenko, and M. van der Schaar. Run-surv: A deep recurrent model for survival analysis. In Artificial Neural Networks and Machine Learning- -CANN 2018: 27th International Conference on Artificial Neural Networks, Rhodes, Greece, October 4-7, 2018, Proceedings, Part III 27, pages 23-32. Springer, 2018.
[22] J. Goeman. penalized: L1 (lasso) and L2 (ridge) penalized estimation in glms and in the Cox model. R package version 09-21 2008, 2008.
[23] J. J. Goeman. LI penalized estimation in the Cox proportional hazards model. Biometrical journal, 52(1):70-84, 2010.
[24] E. Graf, C. Schmoor, W. Sauerbrei, and M. Schumacher. Assessment and comparison of prognostic classification schemes for survival data. Statistics in Medicine, 1817-18):2559- 2545, 1999.
[25] Gui and H. Li. Penalized Cox regression analysis in the high-dimensional and low-sample size settings, with applications to microarray gene expression data. Bioinformatics, 21(13)33001- 3008, 2005.
[26] E. Harrell Jr. K. L. Lee, and D. B. Mark. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine, 15(4):361-387, 1996.
[27] T. Hastie and R. Tibshirani. Generalized additive models: some applications. Journal of the American Statistical Association, $2(398):371-386 1987.
[28] T. Hastie and R. Tibshirani. Generalized additive models for medical research. Staristical methods in medical research, 4(3):187-196, 1995.
[29] H. Hazimeh and R. Mazumder. Fast best subset selection: Coordinate descent and local combinatorial optimization algorithms. Operations Research, 68(5):1517-1537, 2020.
[30] Hill. Employee churn prediction, 2017.
[31] H. Hofmann. Statlog (German Credit Data). UCI Machine Learning Repository, 1994. DOI:
http://dbi.orr//022445//////
[32] T. Hothor, K. Hornik, and A. Zeileis. Unbiased recursive partitioning: A conditional inference framework. Journal of Computational and Graphical Statistics, 15(3):651-674, 2006.
[33] Hothorn, Hornik, and A. Zeileis. ctree: Conditional inference trees. The Comprehensive
R Archive Network, 8, 2015.
[34] C.-J. Hsich, K.-W. Chang. C.-J. Lin, S. S. Keerthi, and S. Sundararajan. A dual coordinate descent method for large-scale linear svm. In Proceedings of the 25th international conference on Machine learning, pages 408 415, 2008
[35] T. Huisman, J. G. van der Linden, and E. Demirović. Optimal survival trees: A dynamic programming approach. arXiv preprint arXiv:2401 04489, 2024.
[36] H. Ishwaran and U. B. Kogalur. Random survival forests for r. R News, 7(2):25-31, 2007.
[37] H. Ishwaran, U. B. Kogalur, E. H. Blackstone, and M. S. Lauer. Random survival forests. Annals of Applied Statistics, 2(3):841-860) 2008
[38] D. Kalbfleisch and R. L. Prentice. The statistical analysis of failure time data. John Wiley Sons, 2011.
[39] J. L. Katzman, U. Shaham, A. Cloninger, J. Bates, T. Jiang, and Y. Kluger. Deepsurv: person- alized treatment recommender system using cox proportional hazards deep neural network. BMC Medical Research Methodology, 18(1):1-12, 2018.
[40] D. G. Kleinbaum and M. Klein. Survival analysis self-learning text. Springer, 1996.

---

[41] J. Ko. Solving the Cox proportional hazards model and its applications. PhD thesis, Master's s thesis. EECS Department. University of California, Berkeley, 2017.
[42] B. Kunzer, M. Berges, and A. Dubrawski. The digital twin landscape at the crossroads of predictive maintenance, machine learning and physics based modeling. arXiv preprint arXiv:2206. 10462, 2022.
[43] M. LeBlanc and J. Crowley. Survival trees by goodness of split. Journal of the American Statistical Association, 88(422):457-467, 1993.
[44] S. Lemeshow, S. May, and D. W. Hosmer Jr. Applied survival analysis: regression modeling of time-to-event data. John Wiley & Sons, 2011.
[45] R. Li and Chambless, Test for additive interaction in proportional hazards models. Annals of epidemiology, 17(3):227-236, 2007.
[46] Y. li, V. Rakesh, and C. Reddy. Project success prediction in crowdfunding environments. In Proceedings of the Ninth ACM International Conference Web Search and Data Mining, WSDM '1, pages 247-256, 02 2016.
[47] J. Liu, Rosen, C. Zhong, and C. Rudin. Okridge: Scalable optimal k-sparse ridge regression. Advances in Neural Information Processing Systems, 36, 2024.
[48] J. Liu, C. Zhong, B. Li, M. Seltzer, and C. Rudin. Fasterrisk: Fast and accurate interpretable risk scores. Advances in Neural Information Processing Systems, 01117760017777 2022.
[49] J. Liu, C. Zhong, M. Seltzer, and C. Rudin. Fast sparse classification for generalized linear and additive models. Proceedings of machine learning research, 151:9304, 2022.
[50] G. Moat and S. Coleman. Survival analysis and predictive maintenance models for non-sensored assets in facilities management. In 2021 IEEE international conference on big data (Big Data), pages 4026-4034. IEEE, 2021.
[51] B. Moufad, P.-A. Bannier, Q. Bertrand, Q. Klopfenstein, and M. Massias. skglm: improving scikit-learn for regularized generalized linear models. Journal of Machine Learning Research, 23:1-5, 2023.
[52] Y. Nesterov et al. Lectures on convex optimization, volume 137. Springer, 2018.
[53] Y. Nesterov and B. T. Polyak. Cubic regularization of newton method and its global performance. Mathematical programming, 108( 108(1): 2006.
[54] M. Y. Park and T. Hastie. I.1 L.-rregularization path algorithm for generalized linear models. Journal of the Royal Statistical Society Series B: Statistical Methodology, 69(4):659-67) 2007.
[55] J. Platt. Sequential minimal optimization: A fast algorithm for training support vector machines. Technical report, Microsoft, April 1998.
[56) S. Pölsterl. scikit-survival: A library for time-to-event analysis built on top of scikit-leam. Journal of Machine Learning Research, 21(212):1-6, 2020.
[57] S. Pölsterl, N. Navab, and A. Katouzian. Fast training of support vector machines for survival analysis. In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML. PKDD 2015, Porto, Portugal, September 7-11, 2015, Proceedings, Part II 15, pages 243-259. Springer, 2015.
[58] T. Popoviciu. Sur les équations algébriques ayant toutes leurs racines réelles. Mathematica, 9(129- (1-1295:2200 1935.
[59] B. D. Ripley and R. M. Ripley. Neural networks as statistical methods in survival analysis. Clinical Applications of Artificial Neural Networks, 237:255, 2001.
[60] M. SÁ Carvalho, Henderson, S. Shimakura, and I. S. C. Sousa. Survival of hemodialysis patients: modeling differences in risk of dialysis centers. International Journal for Quality in Health Care, 1599)189-1196 2003.
[61] R. Sharma, M. Gupta, and G. Kapoor. Some better bounds on the variance with applications.
Journal of Mathematical Inequalities, 4(3):355-33) -33,, 2010.
[62] N. Simon, J. Friedman, T. Hastie, and R. Tibshirani. Regularization paths for cox's proportional hazards model via coordinate descent. Journal of Statistical Software, 39(5):1, 2011.
[63] A. Tarkhan and N. Simon. Bigsurvsgd: Big survival data analysis via stochastic gradient descent. arXiv preprint arXiv:2003.00116, 2020.
12

---

[64] R. Tibshirani. The lasso method for variable selection in the Cox model. Statistics in medicine, 16(4):385-395) 395,, 1997.
[65] V. Van Belle, K. Pelckmans, J. A. Suykens, and S. Van Huffel. Support vector machines for survival analysis. In Proceedings of the third international conference on compulational intelligence in medicine and healthcare (cimed2007), pages 1-8, 2007.
[66] S. Wiseman and A. M. Rush.
learning beam-search optimization. In Proceedings of the 2016 Conference on Empirical Methods in Natural Language Processing, pages 1296-1306, 2016.
[67] S. N. Wood. Generalized additive models: an introduction with R. chapman and hall/CRC, 2017.
[68] B. C.-H. ZHANG. Nearly unbiased variable selection under minimax concave penalty. The Annals of Statistics, 38(2):894-942, 2010.
[69] H. H. Zhang and W. Lu. Adaptive lasso for Cox's proportional hazards model. Biometrika, 94(3):691-700) 2007.
[70] R. Zhang, R. Xin, M. Seltzer, and C. Rudin. Optimal sparse survival trees. arXiv preprint arXiv:2401. 15330, 2024.
[71] Zhu, X. Wang. L. Hu, J. Huang, K. Jiang, Zhang, S. Lin, and J. Zhu. abess: fast best- subset selection library in python and Journal of Machine Learning Research, 23(202):1-7, 2022.
[72] H. Zou and T. Hastie. Regularization and variable selection via the clastic net. Journal the Royal Statistical Society Series B: Statistical Methodology, 67(2)300--320, 2005.
13

---

Appendix to FastSurvival: Hidden Computational Blessings in Training Cox Proportional Hazards Models
Table of Contents
Derivations and Proofs
15
A.1
First, Second, and Third Order Partial Derivatives
15
Partial Derivative y-th Central Moment
20
First and Second Order Partial Derivatives Are Lipschitz-Continoour Continuous
23
A.4
Analytical Solution the Cubie Surrogate Problem
25
A.5
Analytical Solution the regullrrdd Quadratic and Cubic Surrogate Problems
26
B
Related Work
30
Experimental Setup Details
32
Computing Platiorms)
32
Datasets, Baselines, and Licenses
................................
32
Details about Variable Selection Experiments
34
D
Additional Results
35
D.I Optimization on & and f + Eg-rouularized Problems
35
Variable Selection for the CPH Model
43
14

---

A Derivations and Proofs
A.1 First, Second, and Third Order Partial Derivatives
The derivation for the first order partial derivative is shown in Section A.1.1 The derivation for the second order partial derivative is shown in Section A.1.2 derivative is shown in Section A.1.3
The derivation for the third order partial
A.1.1 First Order Partial Derivative
want to show that
(E( 03
Σ
8,
Σ kER;
e
Xx
Σ 8,
ZjER
e
0E(3)
di(n)
0
03
dl(n) Ome
apply chain rule from calculus
mnk
0
tmk 03:
6. logg( Σ t. m/]
# plug in the formula for E(n)
R1
JER,
STR Tk
8
[log(
Σ "
Tmp
03,
jER, # exchange the summation orders; sum over ki, first and sum over i later
8
log(
Tk XKLt
Omp
1=1
jER.
8mk.
evaluate the partial derivative
09
71
0 (log(
0 (n:) XKLl)
8
Σ
e"
# distribute Tiis (-) inside [-]
O1jk1
Omk
jER.
Σ
5,
? 11 Un, 20: 1, Xkit
11= JER ens
1=1
# 0mhi
(log(E
(log(
EJER. e²))) Eyen, en e I 2u. and dmps (1)) 1k
71
Mhs Lyn, Xkil [kI=)XAul)
distribute Xku inside [-]
h1=1 JER,
n 29: Xxyt
Σ
8
U61 I Uin
Σ
1k Xxx1
Σ
EjER,
ens
1=1
=1=1
=1=1
# distribute Σk inside (.)
8
Th Xb,l Σ
l = Xx11)
ybg 2y: J 1
EjER, e
# move 1 1 Zy. the front
E
th Xt,1 Xa)
KYE,, Ejer, e # evaluate the two sunmations over ,y, Le., =i because R2 is set of indices whose time is greater than equal to yi
i the first summation can be simplified
15

---

Σ
8, Σ
1
th Xky1) - Σ
8, X
# distribute Eis 8g inside (.)
MIRR. ZjER. eb
i=1
ΣjER, e Xkii) - -
# move el to the numerator
RI T Xsa)
1=1
8,X
8,
# change notation ky
RE Zjen. eb
A.1.2 Second Order Partial Derivative
We want show that
27
n _
Th
The
8,
X
Σ KER;
Xtt Ejen,
053
Ejen, e%
end
im1
bER
02((B)
(2((T)
033
0
(l(())
OTk2
# Apply chain rule from calculus
0
03
ky=1
0 OTk
ehe
0Tk
Xsa)
Σ
6.Xu
Σ jERR e%
0s:
RER;
i=1
# plug in
at(m))
from the end result for the first order partial derivative above
83t
e
X(a)
Σ
28.X,
TTka RR # distribute Ekz= Oring
EjER, e%
08
amk
0
0 (-)
8762 08,
inside each term inside []; for the first term, also ex-
change the summation orders by summing over i first and k2 later
e
OMka
Xxa)
ΣjER, e%
0
1=1
R.
# second term is 0 because the input to Uika
8
(.) is a constant
e
mmk
Xml
Ejen, emi
03t
1==1
kell.
# move 0 inside the summation of k:
Jmk
EjER, c": (7:)e"h 1 c 2k2 ye,
MTka
X
(EjER, e²))
03,
Ki kEn,
# evaluate the partial derivative using quotient rule from calculus
(EEER c": ee 1==k c 2 yeg veg 2M
Xxa)
(EjEn, eni)
Na v1=1 KER:
8162 #
X
801
(EEER, e"l)) k=
e 1 20, XMIX
(EjER, e²
1=1
kEnl
# move X inside ΣKER.
7 1k=kz
2 c01 1
Way 2014
X Xxa1
ΣjER, en
(EjER, e²/)
divide the fraction into two parts
16

---

e 1 k=ka Xui DEER en
e 20
Nai X
(EjER, e
i=1
N2 #2=1 R.
distribute XuiXuat inside (-)
e limba Xu EjER, en
e Uh2 2u: Xaa Xaa)aa))
(Ejen, e
1=
Ae Ri 1==1
# exchange summation orders; sum over kig and then
e lkmka LjER, en
el yu,2n X
X X
EEjEn, ex) # distribute Ekg=1 inside (+)
KER;
1==1
k==1
e
X
REER, (EER, eno X X X Xrat)]
E
e e
KER; Σ et
for the first term, only k2 k is left because of the indicator kk=kzi i for the second # term, the expression can be simplified because Rí is a set of indices whose time is greater than or equal to ys
X
XuiXuat)
kE Rh E
ZjER, e
(EjER, eni eni)
RER; EER
# distribute ΣKER, inside (.)
X
Σ
Xxl EjeR, e.i
e
(
EjeR, Xust)) e
2
EjER, e
1=1
RER;
kER;
KgER;
# the double sum of the second term can be rewritten product of two sums
X
EjER. eh XN1
XN)
KER; EjeR,
RI
1=1
RER,
# the two terms of the product are exactly the same and can be simplified because k and 22 are independent
A.1.3 Third Order Partial Derivative
We want to show that
3
th Ejen e
Σ9.
82
Σ LER:
X +2
X
Σ
Xxi
0
Ejen, e%
KER
X
Σ MER
Σ EjER, e% XM1 )
ellk
KER, EjER, e/
{²)())
039
0
n
0 (PP((T)
mmks
# Apply chain rule from calculus
DTkka
053
03:
kg
0
T
)²}}
TThs 09
Σ
5
X
XN1)
Omka
Ejen, e
EjER, e
KER.
KER.
plug in
@t(m)
from the end result for the second order partial derivative above
09:
a
Σ
X2
(
e E e X
27 O7ka
mmkss KER.
Σ
EER
03
A=1
KER,
# exchange summation orders; sum over ka and then by moving into the inner summation
a
(-)
00ks inside
8rka
17

---

e
e
27 X
Σ
68
X
X11)
amkss
KER. Σ JER
KER, Σ jER, en
1=1
ka=1
#
Xkal
12
e
Xhi
a Omjies
(E KER;
th EjeR, eni Xks) }
Xxst
JjER edi
hg
RER:
# distribute Dring
a inside [-]
10
0 Oijks
th
a
I
e ΣjER, en Xns 1)
8:
JER:
e)
Orjics
Rg
KER.
KER. Tops
# for the first term, move
()) inside the summation
e
8
X
Otjkcs
EjER
e%
1=1
kg
e
X
e Ejen, eiyy Xxx)
X
jER e%
Drins KER
KER,
# #apply product rule from calculus to the second term
a
e
8:
X2
Ompkss
EJER, ca
1=1
Ny
KER:
Σ KER:
EjER, e% X1))
e
Σ
(
e ESER, en
Xaa X }
amk
NER;
# for the second term, move
(.) inside the summation
12
0 Oikks
58
Xka
EJER, ed
k1=1
KER:
-2
ESER, o% XatXkal
ell
Σ
① Tmk
e ESER. e² Xm1 }}
RER.
KER,
# distribute Xkst inside }
e
elk
8:
1 1k=ks
EEERR e I 2n: X2g
EEER Ed
1=1
kg=1
KER
E
e 1
Xns Xest }} }
Xws)
Ejen, en
EjEn, eT
(EjER, e7/)
2ig 20:
KER.
KER,
9 # Dmss
Eyen, e. Ik=ks
e
e e
1 2m,
ESEn, e"5
(Isen 020)
c els
1 kzka
1 Uky Zw: X2 ka!
ΣjER, e
EjER, e01)2
ky=1
KER. El
e
e elks
XN1
1 =k=
Ejen, em)) Lu 2u Xxs) Xrst
Ejen;
en
JjER, ei
KER;
KER.
KER.
for the last term, distribute
E R. (.) inside [-]
Σ The
e Tha
8
Σ
1 1k=ks
X Xkst
Σ
2 Any Aun zn X XKIX
EjER, ET
(EjER, ens)
kg
KERR
kER,
E
e
e ek
EjER, en XM1
1 1k-ky
Xxl
KER, (EjER, em)) I Vag 2n Xxs) Xho!
Σ
ESER, ens
KER,
KER,
# for the first term, distribute ΣkeR, ()X (xx inside (-)
18

---

8
e
1 kk X
Σ
DJER, e%) 1 20: X X
elt
ZjER,
JER, eb
ka=1 1=1
KER.
MER.
e
XkI
Σ
e
XAs K
DjER, eN
EjER, ken,
en e
R
2
Ejen,
e
en
X11)
EjER, e*
Luny 20 Skal }
kER
# foorr the third term, distribute 2(E KERR Esen, e55 ())))))aan Kugt inside (.)
e
X
1 1k Xkst)
EjER, e
1=1
KER;
k==1
X t
eks
Σ
1 yyn 20 Xkat
Σ EJER.
Σ ΣzER.
KER:
#==1
th kERl
Xxi e
X11)
Σ
1 No Xist)
Σ jen; eN KER
Ejen,
E
k==1
T X
e
20
Zu Xkst]
KER, EjER, en # exchange summation orders; distribute Ekg=1 into each of the four terms inside [-]
EJER, e
k==1
E X
58
Σ
EjER, eni
KER,
X21
Σ
Xkal
Ehg EjER, E9
RER: ZEER, T
kgEli
e
X11
Σ
E EJER, en X
KER, Σ CjER e%
KER,
2 Σ
EN EjER, 0% X1I)
EjER, en Xkst)
e
KER;
RgER;
simplify the summation over K3
E
8:
X
KER EJER, eT
e E JER
e
Xha
X
Σjen, en:
KER:
KER:
t kERl
e!l
XM
X3
Σ jen; eh
Jen, e
KER
e
+
t X11
Ejen, tn Xws)
KER. Σ jen,
KER,
change notation by replacing kg with
n
lk Ejen, eni
8
X
KER:
th
Th EjER. e%
Σ jER, en X2
Σ XM
KER
KER.
E Ejen, eti XM. $1 )
+ 2
# simplify all relevant terms
KER.
19

---

A.2 Partial Derivative of r-th Central Moment
Proof. Recall that the r-th central moment C, is defined as:
th
e
C =
XM Σ
Xxy1
Ljen, en
Σjen, e
RER,
ky EIE
We need to show that
0
(r) = =C++1 - - -0-10
03
0
C+)
03
0
ck
1
Σ
Xyt
Xkal
0
ΣjER,
LAER
JER
ER,
θ
e
Σ
X Σ
Xkal
0,1
ZjER. e
Ejen, en
RER:
KIER;
# move the partial differentail operator 3M
inside Σ R.
e Ejen,
XMI
Xkr!
0
Σjen e
en
KER,
KIER,
T
+
en
XwI E
e
Xhal
Σ
JER, en
0
SJER En
KERR:
# apply the product rule from calculus
XwI Σ
Xk11
0
Ejen, e
Ejen, eNl
KER,
ERR
+
?
r XH
Σ
Xkyl
X
Σ
Xke1 EjER edi
Σ
JERR eti
EjER, eb
03:
KIER,
K.ER.
# apply the chain rule from calculus to the second term L.y" = (r 111y dy
)
(xM-
8
e
Σ KIER
Xka!
0.3
Eyen. eni
EjER,
en
KER.
Σ JER; en
(((xx XMI
r
Σ
EjER, Xkal
e
e
Σ
03,
EjER, e
Tka
Xksl
KIER.
LAER:
# move the partial differential operator in the second term to the front
inside sumn ER; move the negative sign
35
Σ kER:
a
Th
XMI Σ MEE
e Xhy1
0BL
Ejen, e%
Ejen, eT
e
e Ejen, engy X
a e
Xmi
X
Ejen, e
03,
Ejen, e
KER.
ts ERR
LKERR. into the two sunmation terms
# distribute EbER,
ESER, e% XK1l )
Σ RER:
0 e
X Σ
e
0
EjER, 2%
R.
20

---

E
X
Σ
c
Xki
T
Xkyl
0st
EjER,
ens
EjER
en
KIER, EjER, en
ERE
KER,
# move the two values in the second term outside of DKER: because they are independent of
0
e
Σ RER:
XM
Σ
Xkrl
0
E EjER, en
EjER, en MERR,
e
Xki
Cr-1
03, EjER, el
LIE
# siimlify by replacing the last value in the second term with Cr-1 because the central moment definition
Let us focus on the solution for the subproblem a3; 1eth (ESER, e%
0
031
Σ EjER, eds
e
Omk f1
# apply the chain rule from calculus
impk Σ EjER, el
12=1
a amz
OTp B
(c
Σ e
Ejen, e Omic
Ejen, enj
jER
# apply the quotient rule from calculus
2 Ejen, Cnli Iumka
leng22
Xxyt
# calculate the partial derivative
EjER, en
k2=1
E
e
1 kk= Xkyl
EjEn, E )² 1lnng 21, 20, Xkg!
Σ jen, e%
k2=1
k==1
distribbte Eng = (-) Xigt into the two terms inside [-]
2 222
X
Xkvl
evaluate the two summations
EjER, l%
k EjER, enl
c
Xkl Σ
ea Xksl
evaluate the two summations
EjER, e%
EjER,
MER
21

---

Let us now plug this result back into the original problem:
0
(Cr)
03
E
e Σ
Σ
Xki Σ
Xtyt)
0
EjER, en
ela
kE
ME
Xkal
Cr-1
Σ EEER. ens
LiER
# pick up from where we left for the original problem
E ke
e
XHI Σ
e
en Xnss)
(
(XMI-
Σ
e
Xkat
]]
Ejen, e
Ejen,
Σjen,
e
KEER
MEER
r
En, ZjER e%
1
Xkyl Σ
Σ jER, e%
e
Xuzl
Xxsl
Cr-1
kgEll
# #puug in the solution to the subproblem
r+
e
Σ KER,
XxI Σ
XAy1 ens
EJER, EN
EjeR,
MIRR
Σ
Xkal Σ
Xkal
X
Cr-1
LUE EjER, e
EjER, ed
K2ER;
# for the first term, change l2 into ky because both are dummy variables that are independent from each other
r+1
Σ kERl
e
XuI
Σ EIE
XNyt
EJER, ESER tN
Ejen, ens"
X Aurt ] |om
X2
Σ
ZjER e%
ZjEn, jEn, o0
EiE
KsER,
# ffor the second term, move Xbyt inside (.)
r+1
e
Σ keAt
XA1
Σ
X
EjeR, tl
Σjen, e%
kEE
Σ
X2 kil
Σ
2 XhilXnal Cr-1
KIER
JERR
e%
E
R
EjER, e%
for the second term, distribute EnERR EseR,
into the two terms inside (-)
r+17
Σ KER,
e
Xuu Σ
th Xmal
EjER, elt
Ejen en
KIER.
27
2
( Tiq
Σ
X2,
EjeR, en Xhal
C--1
KIER, EjER, e
KIER
in the second term, kn and loz independent dunmy variables, can tumn # double sum of products into a products of sums; we can further simplify this into a square of a because the two terms equal to the same value
C,+1 - T- CC2 Cr-1
# simplify by using the definition of central moment
22

---

A.3 First and Second Order Partial Derivatives Are Supschitz-Continouu
Proof. To show that the first and second order partial derivatives are to show that the second and third order partial derivatives are bounded, respectively.
we need
First order partial derivative is I Lipschitz-conninnuus To show that the first order partial derivative with respect to each coordinate is Lipschitz, we need to show that the second order partial derivative with respect each coordinate is bounded. Recall that the second order partial derivative with respect to the 1-th coordinate can be expressed as:
( a5$
6
XaI)² JER e4
ZJER,
eN
KI
E
RER.
RER:
suffices to show that each term inside the bracket is bounded.
we interpret the expression probabilistically, then the coefficients ESeR, en in front of X2, and Xạ can be thought of as the probability a particular distribution because all terms are greater than or equal to 0 and sum up to 1.
For notational convenience, let us use a to denote the probability of this specific distribution with Then we can rewrite each term inside [+] as:
ag
Σ Ejen, 5
Σ RER,
e
X
Σ kER.
E e Ejen.
Xm1 X=)²=
Σ x E Xx))
DjER, e
kER
kER.
The right-hand side is nothing but the variance of {XKINEER, with respect to the distribution JKER. Since the variance is always non-negati non-negative, have
Σ X
E XuI) 20
RER,
Let us now denote = minkeR R, X&t as the minimum of this given set, b b: MAXKER, Xur as the maximum of this given set, Z as a random variable with values restricted to [a, bj. We are going to show that
Σ Eax
Σ xxI) max -EE()
(23)
RER:
R
We achieve this through two steps.
First, suppose the random variable can only take finite |R number of values {Z1. Zz...., .... (R. with probability (P1, P2, PIR,1). where |R is the cardinality of the set R, Then, have
Σ X - C xxxx
max [axX[E,|Z? (E,!2)))
(24)
Z.P
kEn;
kEN;
where the expectation is taken with respect to the distribution p. The above inequality holds because the left-hand side is a specific instance of the expression inside [t]. so taking max(-) produces the inequality above.
Next, if we drop the assumption that the random variable Z can only take [&, number of values, then the maximum variance can achieve is smaller than before. Mathematically. this means that
max E | 7.p
(
max E z (B[Z))
(25)
Combining Inequality (24) and Inequality 25 we arrive at Inequality (3).. Lastly, note that the
Popoviciu Inequality [58] tells us that for bounded random variable restricted to [a, b). the maximum variance it can achieve (a-b)² i.e.
(b-a)² a)²
max E E
(
This allows us to conclude that
bb a)²
(maxken, XAt - minken, Xaa)
Σ X2i
Σ
Xxl
<
ke
WE
R.
23

---

Plugging in this inequality into the expression of the second order partial derivative, we can show that the second order partial derivative bounded:
)E((3) 0.9
e
0 <
8
X
Ejen, e Xx,1)
ZJER, JER
e
i=1
E
(max KER;
Xxt
min XRa)² KER;
Second order partial derivative is Lipschitz-continuuss This proof is similar to the first part above. We need show that the third order partial derivative with respect to each coordinate bounded. Recall that the third order partial derivative can be expressed as:
th Xxa)²
8
X + 2(
Σ
Σ JER, e
E JER e
1=1
KER,
KER;
e
X
Zjen, e X
Ljen,
t
KER.
KER.
As we have done in the first part, if we to denote the probability of this specific distribution with ea
Esen, e we can rewrite each term inside [=]
Σ
t
X +
Σ
T
XM1
-3
e
)x ken,
X
Xal
EjER, eT
EjeR, e%
EjeR, E
Σ
kEn;
kER;
kER,
Σ X + 2( KER
ag Xx)3
3
Σ a X
Σ axXMa
KER.
KER.
kER.
Similarly, as in the first part, let us now denote a minsER, X& as the minimum of this given set, b ma Xa as the maximum of this given set, Z random variable with values restricted to (a, b]. Using the exact same logic, we have
Σ KER,
ag X + 2(
0g. XNI)
3( Σ
Og X
Σ ag Xx)I
R
kE R.
KE
max 3| + 2E 2 3
3E [E
-mxx [[[[ -E
The expression inside |-| the right-hand side is known as the third central moment (skewedness) in statistics. Fortunately, we can derive an explicit formula for the maximum of the absolute third central moment of bounded variable.
According to [61], we have the following inequality involving the second and third central moment:
[E[ - 3|
E - = +
S (b a)²
2E - 2
From this, we can derive an upper bound on the third central moment:
[Z - E| < (Z 1|
(6 a)² E[(Z - E
For notational convenience, let denote V := E[(Z ] then the right-hand side above can be expressed as a function V:
f(V) = 2V V 7 a)² == Vbb-a)222 - -VV3
Because V is the variance, V E [0, (b )) Additionally, because /i is monotonically increasing, the maximum is achieved either at the points where the first order derivative of (b a)² 4

---

with respect to V is 0 or at the boundary, 0 and {bb - -a)?. Let calculate the points where the first order derivative
d dV [b a)²v² 4V 2(b - =
(b aa²
V = 0
To obtain the maximum value achievable, we calculate the values at points V = 0, V (b-a)² V = (b - a²² and pick the maximum value afterwards:
and
6
f(0) 2 01 1 a)² -
=0
1 1 a)9) 2 (b a)²
1
1
1
V a)²
(b (b-a)- 0)²
al
6
6v3
f1 a)²) 2 x
1 (6 a)² a)²
1
(b a)² 0
4
Therefore, the maximum value achievable is
5/3 al3.
We now show that upper bound on the absolute value the third central moment is actually tight by providing with a concrete example. For a random variable Z, let IP|Z = aj 1. PIz b] = 2. and
P|Z ا 1.2 We can verify that E[(Z B[Z])*
o b al3, thus proving that this upper
bound is indeed tight.
This helps us to arrive at the following inequality:
1 max Xai min Xaa KER;
Σ X kt
2(
Σ ag Xm)²
ag
Σ
6v3
KER,
EER,
KER,
kER;
Therefore, for the upper bound of the third order partial derivative, we have the following explicit formula:
8. max XMI
min XRsl 3 RERt
033
6v3
kERl,
A.4 Analytical Solution to the Cubic Surrogate Problem
Let f() be convex function whose first, second, and third derivatives all exist. Let h (Ax))) = f(x)) + A*) (x) 1 Ar/ be the cubic surrogate function (53) of [(x,, where
L3 > 0 is the Lipschitz-constant of the second derivative fx()). Then, the minimum of this surrogate
function is achieved at the following point:
Ai argmin hn (Ax) sgn(f((x)) f"(x) (T)( + 2Ls (())
(26)
Ar
L3
Proof. We discuss three cases: f'(x) > 0. (z) (z) < 0, and f'(x) =0
Case f'(x) 0.
If (x) > 0, then Ai < 0. For the sake of contradiction, suppose A. > 0, which means h (Ax)
achieves its minimum some point with Ar 0. However, we arrive a contradiction because
hr (0) fx) < f(x) (r) + r" Ln h (Ar), for Ax 0.
Therefore, the minimum is achieved either at Ar = 0 or Ar < O. However, since Ar and |rr| grow slower than |Ar when Ar. is close to O, there exists some Ar 0 such that f'(x) Ar f Az) Az < O. Thus, (0) cannot be the minimum value, and are left with the
25

---

minimum value achieved at some A < <. Ar < = - f())+ (x)))== ( Note that the second order derivative of o is greater than equal to 0 since
(()) Lg 2 0.
dA 2 O because [(x) is convex, 0 S 1 (x)) S Lag and Ax < 0
Therefore, hz (Az) is convex function with respect to Az when Az < 0, and its minimum value is achieved when the first order derivative is O. When the derivative with respect to A is O, we have
[Ax)) =0 0 Ar (x) (f (2)) (z)
f(x)+
2
L3
(x())- ((<))V(("()) ((M)))(() satisting the condition Ar 0.
Since f(x 0, only one root
Lg
Thus, when (x) O, we have
As (* - (z)) (z)
L3
Case 2 (x) < 0.
If fx)) 0, we have A& > using the same logic as above. If Ar 0, (Ar) f(x)+ f(())x+ + (z) + (Lx)) Note that the second order derivative of (z(Az) is also
greater than or equal to 0 since
(Ar) = (x) La 2 O.
dA
# (f) because (xx is convex, L3 > O, and Ar
Therefore, h2 (Ax) is convex function with respect to Az when Az > O, and its minimum value is achieved when the first derivative When the derivative with respect to Ar 0, we have
f(xx) + (A))) - (x) = (z) (z)
L3
-x'(x) ((')))) ((("))))22223"1 (L)
Since f'(x) < 0, only one root
satisfing the condition Ar 0.
Lg
Thus, when fxx 0, we have
A -f" + V (+) (z)
L
Case 3 fx(()) 0.
When f'(x) the minimum value of b2 hr (Az) is achieved Ax 0.
The explicit formulas for the three cases above can be unified into one succinct formula below:
A sgn(f((z)) ('(x) (ff (x))² (f)
L3
A.5 Analytical Solution to the & f--reuuarized Quadratic and Cubic Surrogate Problems
{ regggularde quadratic surrogate problem surrogate problem:
We have the following & l-rreularrzed quadratic
Ab argmin axx bA + +11e Ar).
Az
The solution for the above problem
a \1)/b if be a -
A6
a + X1)/b if bc- > d
otherwise
26

---

Proof. Since the function aAx + =bA²² + A1 Arl is convex, the condition for this function to achieve the minimum value is for its differential include 0. The differential of this function is:
a dAr A1
if Ar > -c
Os. (ade + joso +11c+ + Sel)
a a bAr A1
if Ar < C
a a b&r 11 X1]
if
Ar = -C
1. For the first condition, if the differential contains 0, we have
Ar +
a + bAx + A1 =
b
However, because we require Ar > -c, we have
- A1
2. For the second condition, the differential contains 0, we have
a - A
+ dAr - A1 0 Ar=
b
However, because we require Ar < -c, C, we have
0 A1 -c bc < ー
3. For the third condition, if the differential contains O, we have
0 E a + bår + [->1, A1] = a + b&x A1 SoSa+bArth
However, because we require Ar cc we have
-λ1 s bc a \1
é--reuularized cubic surrogate problem problem:
We have the following & f--reullarized cubic surrogate
Aa argmin aAz + 2 =bd0a ccl0 + +1ild Ar/.
(27)
Az
The solution to the above problem is
sg + / 2e(sgn(d)) + +1)
/c sgn(d)a + A1 $ 0
sgn
b v 82 2e(sgn (d)a X1)
/e sgn(d)(a bd) ed > A1
A#=
sgn(d)
b+ /B2
(d)a + \1)
/c
if sgn(d)(a bul) 1 cdf < ->1
-d
otherwise
Proof. Like the first part, since the function aAr + 1.b22² + f + Aj + Ar/ convex, the condition for this function to achieve the minimum value is for its differential to include O. We discuss the differential of this function in two cases: 2 0 and d < 0.
When d 2 0, the differential of this function is:
ds
c// + AjId + Ax| aAr + 7 bax
b
2 + A
if Ar > 0 if Ar = 0
c Ax? Ao + A1
coA + A
if
d < Ax <
a AAT
ccA + 1->1;>i
if Ar = -d
bAz
1cA - A1
if Ar < -d
We discuss these 5 cases one by one.
27

---

1. For the first condition, if the differential contains 0, we have
Ae -b = /b² 2c(a + X1)
a b c + A1
However, because we require Ar > 0, we have
\1 (a + bA + -
This means that can only have one root because the other root violates Ar 0:
A b + 12 2c(a + X1)
2. For the second condition, if the differential contains O, we have
0Ea+ bba + 2 A A+ + X
aa+ bAs c + A1 S 0 < a 6 + c X
However, because we require Ar = 0, we have
a + \ = 0
3. For the third condition, if the differential contains 0, we have
DAt cAi + A 0 A#= + 02 + 2c(a + A1)
However, because we require
-d Ax < 0, we have
A = a bAt
c > >- >
2
This means that we can only have one root because the other root violates the condition d < Ax 0:
A b D2 2c(a λ1)
C
Moreover, since the root is between -d and 0, and the coefficients in front of Az. -/c. negative, have a bd tedl A 0. 4. For the fourth condition, if the differential contains 0, we have
0 + AAi 2
+1-111
+
1
a bAs However, because we require Ar = dd, we have
cA2² 505 sa bAd
- A1
2
a bd
2 od A
bd od A1 la bd
1 off <
For the fifth condition, the differential contains O, we have
b c >i + A -b + 62 + 2c(a - A1)
However, because we require Ar < d, we have
\1 a bAg 1
cA2² <a
This means that we can only have one root because the other root violates the condition Ax < dd < 0:
Ai b 1 + 2c(a A1)
e
Moreover, because the left root is less than -d d and the coefficient in front Az?, 2c is negative, we have bd ced > 0.
1
28

---

When d < 0, the differential of this function is
a + bar
+
7c// A1 + Arl)
1
b
cA + X1
if Ar > -d if Ar -d if 0 Az <d if Ar=
bAr Ccr2 [-11;
bAT + cAx2 \
+ dA + +2-+t---000 c/-Ax² "r,1 1
X1
+ bAT ccA - X1
if Ar <<<0 0
Similar to the previous part when d 2 0, discuss the 5 cases one by one but omit the details because the logic and the reasoning process are exactly the same:
For the first condition, if the differential contains 0, we have
A b T2 2c(a A1) and bd + Jol
+ 11 <0
2. For the second condition, if the differential contains 0, we have
A -d and |a bd cd S A1
For the third condition, if the differential contains 0, we have
A b + 82 20(a X1)
4. For the fourth condition, if the differential contains 0, we have
Ai = 0 and Åι =
5. For the fifth condition, if the differential contains O, we have
1 2c(a - A1) and a 2
A
c
We can combine the two situations where 2 and d < 0 and obtain a unified formula:
sgn(d) + / 2c(sgn(d)a /e
sgn(d)a + A₁ <0 0
sgn(d)a + a1 = 0
As
sgn(d) b b+ v + +e(eesgn())de X1) /e
sgn(d) - bd) Ffd > A1
sgn(d)
b / 62 20(sgn (d)agn(dda>>)) - + 21) /c
sgn(d)(a bd) of < -
-d
otherwise
The first and second equations above can be further unified into just one equation sgn(d) (-b
0 2c(sgn + A1)) if sgn(d)a + A S 0. we finally have
sgn(d) [b + /B 20 (e((sg((d))) X1)
/c sgn(d)a + \1 5
sgn(d)
b+ /b2 +22 d)a X1)
/c if sgn(d)(a bul) ced A
Ai
sgn(d)
b+ / 20 (da 21) /e if sgn(d)(a bd) cfl < -
-d
otherwise
29

---

B Related Work
Optimization for CPH One way to train the CPH model is through gradient descent 63 However, because the CPH loss is complex, it is difficult to pick the right step sizes, so gradient descent tends to be slow when we want to obtain solutions with high precision. To speed up the training process, people have used the Newton method [25] 22 23]. The drawback of this approach is that is computationally intensive to calculate the full Hessian matrix. Moreover, the vanilla Newton method cannot be applied to solve the f--regularized problem. To alleviate this problem, quasi Newton [62] and proximal Newton [51] methods have been proposed. However, as we have shown in our experiments, these Newton methods have the flaw of training loss blow up due to vanishing second order derivatives. A generic binary search method [4]] has also been proposed, but the algorithm has been shown to be slower than the quasi Newton method. In contrast to all these approaches, our method is computationally efficient, can easily handle different regularizers, and guarantees that the loss decreases monotonically.
Modern First and Second-order Optimization Methods Our work is greatly inspired by the modern developments for convex optimization [52]. but the general principle cannot be rigidly applied. For first order methods, as we mentioned about gradient descent in the last paragraph, it is difficult choose the right stepsize for fast convergence. Instead of performing gradient descent, our method performs coordinate descent, which has been shown to be effective in training other statistical We give an explicit formula, by leveraging the Popoviciu's 's inequality on variances 158 to calculate the Lipschitz constant at each coordinate. We also design a second order method (still under the coordinate descent framework) based on the cubic-regulariization of the Newton method [53]. To calculate the Lipschizt constant, we make the connections to the third central moment and Shatia-Davis's inequality (61). Moreover, we are able to exploit the mathematical structures of the CPH model to compute the second order partial derivatives at the computational complexity of On),, making the evaluation per iteration of our second order method as fast as that of our first order method.
models 9
18
19
20
34 55)
Variable Selection and Interpretability for the CPH Model If we can find sparse solutions whose predictive performances are as good as dense solutions, we can better interpret which features play important roles. A popular way to select important variables for the CPH model adding an & regularization term (64), commonly known as the LASSO method. We can also use the FlasticNet method, which adds an & + 82 regularization term 622. Another way to apply the Adaptive LASSO [69], which repeatedly use the absolute values of coefficients obtained from the previous iteration as weights of & regularizations for parameters in the current iteration. These above approaches all use convex regularizers and have difficulty obtaining high-quality solutions when the support size is small. The reason is that these convex regularizers penalize the magnitude of the coefficients while promoting sparsity. To avoid this issue, recently, solving the fg-constrained problem [49] [29] 4 has attracted lots of attention and shown to produce much sparser models without losing accuracy. For the fo-oonttrddeed CPH problem, ABESS [71] has proposed to hybrid method of greedy selection and feature swapping to solve the problem heuristically. However, as have shown in our experiments, this method cannot handle highly correlated features. Our method also solves the fo-constrained CPH problem but uses the heam search framework [66] 48 47). Although the beam search framework has already existed, this frameowork cannot be applied the CPH model without our coordinate descent algorithm. We need to use coordinate descent for support expansion as well as coefficient finetuning, in which other Newton-type methods all have issues with losses potentially blowing up.
Other Model Classes for Survival Analysis In addition the CPH model, there are some other model classes that can be applied to analyze time-to-event data. One model class is the survival tree models 70 3 35 Survival trees have the advantage of capturing non-linear interactions between features. However, when sparse models are desired, the accuracy of sparse trees is compromised by the fact that all samples in the same leaf node share the same predictions. One way to overcome this issue is to construct ensembles of trees using random forest or boosting techniques [6) 32 33
Another model class for survival analysis is based on neural networks 39 6
59 21. However, for
all these other model classes mentioned, they are not very interpretable due to large parameter space. The CPH model, which is the focus of our work, provides both interpretability and good accuracy. For applications involving high stakes decisions, it is desirable to produce models that are as sparse
30

---

as possible without losing accuracy. In this work and especially the variable selection experiments, we push the limit of pparsity-accuracy tradeoff curve for this model class.
31

---

C Experimental Setup Details
C.1 Computing Platforms
All experiments were run on the Intel(R) Xeon(R) CPU E5-2680 v3 Processor, 2.50GHz. We set the memory limit to be 100GB.
C.2 Datasets, Baselines, and Licenses
We have a summary datasets for experiments in Table
Dataset
Samples
Origin Features
Encoded Binary Features
Flchain
7874
39
333
Kickstarterl
4175
54
2144
Dialysis
6805
7
207
Employee
14999
17
272
Synthetic HinthtticfighCCoorHiiiiiii
1200
1200
N/A
SyutheticHighhCrrHgghbD
900
900
N/A
Syytheeic
600
600
N/A
Table 1: Datasets Summary.
Synthetic Data Generation Process the following process, similar to [71]:
The synthetic data used in the paper is generated according to
Firstly, from a Gaussian distribution (O,, Σ) where the first entry the mean and the second entry is the covariance matrix with size p x p, we sample features:
a; ~ MO, E).
(28)
The covariance matrix is defined as Ejt = اااااا where E (0, 1] is the correlation parameter. When P is large, the features x; are more correlated. We create k-sparse coefficient vector 3 RP. The entries of 3 are either or 0. If (p/kk//)) = 0, then 33 = 1; otherwise, 33 = 0.
Secondly, generate the death time f; according to the following equation:
= (- log V/ c 3
[olog
(29)
where V; ~ (O, 1) (samples are drawn from uniform distribution on the interval [0, 1]) and is hyperparameter. our experiments, we set s = 0.1.
Lastly, we generate the censoring time, the censoring indicator, and change the death time to observation time. We sample the censoring time from uniform distribution: ~ (O, 1). the
death time is bigger than the censoring time, we have the indicator equal to 1; otherwise, we have the indicator equal O. Specifically, we have:
6, = L4>c
(30)
Afterwards, we change the death time to observation time, taking into consideration of censoring:
t; min(t;, Ct).
(31)
We form a triplet (æ,, 4 6.) and return this triplet as one sample.
Real-world survival data:
Use of nonclonal serum immunoglobulin free light chains to predict overall survival in the general population [00. The event death.
Flchain:
* Kickstarterl:
Data from popular crowdfunding platform, used to predict project suc-
cess (46). We used the version from https:: :tttps://dmkkd es. wtt.eed//rrjjecss/urriii
data/
32

---

Data from a survival study of dialysis patients, which aims to assess quality of renal replacement therapy dialysis centers in Rio de Janeiro, Brazil [60]
Dialysis:
EmployeeAttrition: The task of predicting when an IBM employee will quit. The event is IBM employee's leaving (30]
Licenses
We list the licenses of the software packages used in this paper:
Abess: The license of this package is GPL-3.
skgim:
The license of this package is BSD-3.
Scikit-survival (SkSurv): The license of this package is GPL-3.
Flchain: We use the dataset from the Scikit-Survival [56] package. The GitHub
link this dataset https://gt cwmp
sksurv/datassets/daaa. The license of this package is GPL-3.
Kickstarter1: We the dataset from the Virginia Tech. The link to this dataset is There no license associated with this dataset. This means cannot modify any part of the dataset, which we have obeyed while doing experiments on this dataset.
hattps://mkk.. cs vt www.eeduuoooooo pttppptssssrrriaaalllla
Dialysis: We use the dataset from the SurvSet 12 package. The GitHub link
this dataset is httppp:://aa..gttttttueeeeeeooos https///seaaaaaaaoooooooo////a/a httpee///daaaaggnoouuuuttt///bbe
csv The license of this package is GPL-3.
EmployeeAttritton We the dataset from the PySurvival 17
package. The GitHub
link to this dataset is http:///wwb..coooggaar////rrrraa
atw.riiinnttttiiiioetooo cav The license of this package is
Apache-2.
Baselines
We compared our method against various survival models:
Abess:
Adaptive Best-Subset Selection (ABESS) algorithm [7] for Cox proportional hazards model. We used the Cox model in abess python package Version 0.4.6.
SkglmALassooox: Cox model with the adaptive Lasso regularization 69
We used the
implementation from skglm [ [1
SksuryCoxnet: Cox's proportional hazard's 's model with elastic net penalty 622 We used the implementation from Scikit-survival (SkSurv): seikit-survival version-O.20.0 (https
readthedocs. 10/ /ttp///////l//iee......e html).
SksuryTree: A greedy decision tree model using log-rank splitting rule
43 We used
implementation from sksurv.
SksurvRSF: Random survival forest 37 algorithm. We used the implementation from Sksurv.
SKSurvGBST: Gradient-boosted Cox proportional hazards loss with regression trees base learner. In each stage, a regression tree is on the negative gradient of the loss function. We used the implementation from Sksurv.
SksurvNaiveSVMM: Naive version of linear Survival Support Vector Machine the implementation from Sksurv.
65
We used
SksurvFastSVe Efficient Training of linear Survival Support Vector Machine used the implementation from Sksurv.
57
We
Evaluation Metrics
The full name of this metric score is Harrell's Concordance Indices 26 is used to evaluate the discrimination ability of a survival model. It assesses how well the model ranks observations based on their predicted risk of experiencing an event (e.g., death, discase recurrence) over time. The higher the CIndex score, the better the model.
CIndex:
33

---

2. IBS: The Integrated Brier score was proposed by [24] to evaluate survival models across all possible time threshold. The IBS score takes the Brier score a step further by integrating across all possible time points within the follow-up period of interest. This provides single score summarizing the model's performance over the entire time range. The lower the IBS score, the better the model.
3. Fi-score, Precision, Recall: Suppose the true coefficients are 3r and the estimated coefficients are B. Then the precision score can be calculated as P = |supp(3*) U (ppp(3)/\supp($})) where supp(-) extracts the support (indices whose coefficients are nonzero) of the input vector. The recall score can be calculated as R (supp((*) supp(3)///sppp($*) We calculate the Fl score as Fl = 2PR/(P + Rt..
C.3 Details about Variable Selection Experiments
We ran 5-fold cross-validatio (random seed 0) on the following datasets: Dialysis, Flchain, Kickstarterl, Emplloee SyythetieHighCCorHiggiiiii SyntheticHigh HHighDim2, HighheeiciihCCoHHgghhh continuous features into binarized features, by considering 1000 quantiles for each continuous column. For each dataset, we ran algorithms with different configurations and evaluated fitted models with metrics described in Appendix C.2
Collection and Setup:
In order to create highly correlated features, we encoded
Abess:
We this algorithm with 30 different configurations: support size, k ranging from to 30, forcing the number of non-zero coefficients in the Cox model to be exact k. We set primary_modes fit__ iit_mar__ae jter to be 20, approrimate Newton be False. All other parameters were set to the default.
SksurvCoxnet: We ran this algorithm with 30 different configurations: support size, k. ranging from 1 to 30, forcing the number of non-zero coefficients in the Cox model to be exact k. We set [1_ [1_ratio be 1.0, alpha_mite mai__mnnnnmom r_tiio be 0.01. All other parameters were set to the default.
LaasmmALsssoComm We ran this algorithm with different I.1 regularization penalty param- eters (alpha): 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100. All other parameters were set to the default.
We ran this algorithm with 8 different configurations: max depth limit, d, ranging from to 9, and corresponding maximum leaf limit 2d. The random state was set to 2024 and all other parameters were set to the default.
SksurvTree:
SksurvRSF, SKSu We ran this algorithm with 5 configurations: max depth limit, d, ranging from 2 to 9, and 5 different total numbers of estimators (10, 50, 100, 500, 100). The random state was set to 2024 and all other parameters were set to the default.
MisurrFastSVM. We ran this algorithm with 9 different &g regulariza- tion penalty parameters (alpha): 0.01, 0.05, 0.1, 0.5, 5, 10, 50, 100. All other parameters were set to the default.
SksurvNaiveSVM.
SeSurrCCuPHBeemSeaar (our method): We ran this algorithm with 30 different config- urations: support size, ranging from 1 to 30, forcing the number of non-zero coefficients in the Cox model to be exact k.
Recording experimental results: For each method with specific configuration, we have a set of up to fitted models on each dataset. Some metrics may be unavailable:
Precision, recall, and fl-score are not available on real-world data as we do not know the true coeffcients.
The losses on the training and testing folds of cox models are not applied to non-Cox models.
SksurvNaiveSVM and SksurvFast VM were not able to provide IBS and AUC.
Some methods' training time exceeded our 3-hour time limit (We noticed that sksurv- NaiveSVM often timed out).
For Cox and SVM models, we recorded the number of non-zero coefficients as the support size. For tree based models we recorded the number of nodes as the support size. We plotted the standard deviation of support size and various metric scores corresponding error bars.
34

---

D Additional Results
D.1
Optimization on & and ( + fz-regularized Problems
D.1.1 Results on Flchain
Loss Iters
vs Time
Dataset: Fichain
Dataset: Fichain
10
losses
10a
10a
10
15
20
0
iters
elapsed times (sec)
NewtonExact
NewtonProx
NewtonQQasi
QuadraticSurrogate (Ours)
CubicSurrogate (Ours)
Figure 5: Optimization on the Flchain dataset with A 0 and X2 1.0. The baselines (exact Newton, quasi Newton, and proximal Newton) all have the losses blow up. contrast, our methods based on the quadratic and cubic surrogate functions have the losses monotonically decreasing.
Loss Iters Dataset: Fichain
vs Time Dataset: Fichain
600000
losses
300000
100000
4
10
12
14
0
5
10
15
20
iters
elapsed times (sec)
NewtonExact
NewtonProx
NewtonQuasi
Cuedraticsurropete (Ours)
Cubicsurogaate (Durs)
Figure 6: Optimization on the Flchain dataset with A = 0 and A₂ = 5.0. The baseline, exact Newton, has the losses blow up despite stronger & regularization. The other baselines, quasi Newton and proximal Newton, do not have this issue when &g increases but are significantly slower than our methods.
35

---

Iters
vs Time
Dataset: Fichain
Dataset: Flchain
10"
16a
losses
10
20
0
1
2
3
iters
elapsed times (sec)
NewtonProx
NewtonQuesi
QuadraticSurogatte (Ours)
CubicSurregate (Ours)
Figure 7: Optimization on the Flchain dataset with A 1.0 and X2 = 1.0. The exact Newton method can applied solve the --rreuulaiieed problems, and proximal Newton. These two baselines both have the losses blow up when the l2 regularization weak. In contrast, our methods based the quadratic and cubic surrogate functions have losses that monotonically decrease.
only compare with quasi Newton
Iters Dataset: Fichain
vs Time Dataset: Fichain
15000
15000
14500
14500
losses
13500
13500
13000
13000
12500
12500
12000
0
2
10
12
14
0
2
8 10
12
iters
elapsed_
pooosed_tmmess times (sec)
NewstnnProx
NewtonQuasi
QuadraticSunogate (Ours)
CublcSurrogate (Ours)
Figure 8: Optimization on the Flchain dataset with A 1.0 and X2 = 5.0. The exact Newton method can be applied solve the & ferreuularized problems, only compare with quasi Newton and proximal Newton. Stronger & regularization helps these two baselines avoid the losses going into infinity. However, these two baselines still significantly slower than methods.
36

---

D.1.2 Results on Employee Attrition
Loss Iters
Loss Time
Dataset: EmployeeAttrition
Dataset: EmployeeAtrrition
24500
24500
24000
24000
23500
23500
losses
22500-
22500
22000
22000
21500
21500
10
12
0.0
2.5
5.0
10.0
12.5
15.0
17.5
iters
elapsed_times Limeedtttees_tteeeee (sec)
NewtonExact
NewtonProx
NewtonQuasi
QuadraticSuroogate (Ours)
CubicSurogate (Ours)
Figure 9: Optimization on the Employee Attrition dataset with è and A2 = 1.0. Although our methods make less progress toward the minimum loss per iteration (left plot), faster than other methods terms elapsed time (wall clock) due to cheap evaluation cost per iteration. For ease figure reading, only give partial plot with few iterations. When the number of iterations is large, our methods achieve better losses than the baseline methods.
significantly
Loss Iters Dataset: EmployeeAttrition
Loss Time Dataset: EmployeeAtrrition
24300
24000
23500
23500
losses
23000
23000
22500
22500
22000
21500
21500
10
12
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
iters
elapsed times (sec)
NewttoEExact
NerwtonProx
NewtenQuussi
Quedraticsurrogate (Ours)
Cubiesurogete (Ours)
Figure 10: Optimization on the Employee Attrition dataset with A and Xı 5.0. Although our methods make less progress toward the minimum loss per iteration (left plot), we are significantly faster than other methods in terms of elapsed time (wall clock) due cheap evaluation cost per iteration. For ease of figure reading, we only give partial plot with few iterations. When the number of iterations large, our methods achieve better losses than the baseline methods.
37

---

Iters Dataset: EmployeeAttrition
vs Time
Dataset: EmployeeAtrrition
24500
24500
23500
23500
losses
23000
23000
22500
22500
22000
21500
21500
10
12
0
i
6
iters
elapsed times (sec)
NewtonProx
NewtonQuesi
QuatraticSurogatte (Ours)
CubieSurrogate (Ours)
Optimization on the Employee Attrition dataset with À = 1.0 and À2 = 1.0. The exact Newton method can be applied to solve the f-regularized problems, so we only compare with quasi Newton and proximal Newton. Although our methods make less progress toward the minimum loss per iteration (left plot), clock) due to cheap evaluation cost per iteration. For ease of figure reading, we only give a partial plot with a few iterations. When the number of iterations is large, our methods achieve better losses than the baseline methods.
Figure 11:
significantly faster than other methods in terms of elapsed time (wall
Loss VS Iters Dataset: EmployeeAttrition
Loss vs Time Dataset: EmployeeAttrition
24500
24500
24300
24000
23500
23500
losses
28000
23000
22500
22500
21500
21500
iters
elapsed times (sec)
NewtonProox
NewtonQuesi
QuadraticSurrogate
CubicSurrogate (Ours)
Figure 12: Optimization the Employee Attrition dataset with A 1.0 and A₂ 5.0. The exact Newton method can be applied to solve the & reegullareeed problems, so we only compare with quasi Newton and proximal Newton. Although our methods make less progress toward the minimum loss per iteration (left plot), we are significantly faster than other methods in terms of clapsed time (wall clock) due to cheap evaluation per iteration. For ease of figure reading, only give partial plot with a few iterations. When the number of iterations is large, our methods achieve better losses than the baseline methods.
38

---

D.1.3 Results on Kickstarterl
Loss VS Iters
Loss Vs Time
Dataset: kickstarter1
Dataset: kickstarter1
losses
10a
10a
10
15
40
60
80
Iters
elapsed times (sec)
NewtonExact
Newtonnrox
NewtonQuasi
QuedraticSurrogate (Ours)
CubicSuriogate (Ours)
Figure 13: Optimization on the Kickstarterl dataset with A = 0 and X2 1.0. The baselines (exact Newton, quasi Newton, and proximal Newton) all have losses that blow up. In contrast, our methods based on the quadratic and cubic surrogate functions have monotonically decreasing losses.
Loss Iters Dataset: kickstarter1
vs Time Dataset: kickstarter1
103
10a
losses
10f
10
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
0
20
40
60
80
100
120
iters
elapsed times (sec)
NewtonExact
NewtonProx
NewtonQuasi
QuedraticSurrogate (Ours)
CubicSurrgaate (Durs)
Figure 14: Optimization on the Kickstarterl dataset with A = 0 and X2 5.0. The baselines (exact Newton, quasi Newton, and proximal Newton) all have the losses blow up. In contrast, methods based on the quadratic and cubic surrogate functions have monotonically decreasing losses.
39

---

Iters Dataset: kickstarter1
vs Time
Dataset: kickstarter1
6-107
6+33
44127
losses
1+197
2412
10
0.0
2.5
5.0
7.5
10.0 iters
12.5
15.0
17.5
20
40
60
80
elapsed times (sec)
NewtonProx
NewtonQuesi
QuatraticSurogatee (Ours)
CubieSurrogate (Ours)
Figure 15: Optimization on the Kickstarterl dataset with A = 1.0 and A = 1.0. The exact Newton method can applied solve the & rerreuulri problems, and proximal Newton. These two baselines both have the losses blow up when both the & and l2 regularizations weak. In contrast, our methods based the quadratic and cubic surrogate functions have monotonically decreasing losses.
only compare with quasi Newton
Iters Dataset: kickstarter1
Vs Time Dataset: kickstarter1
300000
250000
250000
200000
200000
losses
150000
150000
50000
50000
0
0
2
4
6
8 iters
10
12
14
16
0
20
40
60
80
100
120
elapsed Limees (sec)
NewstonProx
NewtonQuasi
QuadraticSunogate (Ours)
CublcSurrogate (Ours)
Figure 16: Optimization on the Kickstarterl dataset with A₁ = 1.0 and X = 5.0. The exact Newton only compare with quasi Newton
method can be applied to solve the l rerelullrrd problems, and proximal Newton. These two baselines both have the losses blow up even when we have stronger & regularization. In contrast, methods based on the quadratic and cubic surrogate functions have monotonically decreasing losses.
40

---

D.1.4 Results on Dialysis
Loss Iters
Loss Time
Dataset: Dialysis
Dataset: Dialysis
10300
10300
10200
10200
losses
10100
10100
9900
9900
9600
2
5
6
iters
elapsed Liteedtttese_teee (sec)
NewtonExact
NewtonProx
NewtonQuasi
QuadraticSurrogate (Ours)
CubicSurrogat (Ours)
Figure 17: Optimization on the Dialysis dataset with X and A2 = 1.0. Although our methods make less progress toward the minimum loss per iteration (left plot), we significantly faster than other methods in terms of elapsed time (wall clock) due to cheap evaluation cost per iteration. For ease of figure reading, only give partial plot with few iterations. When the number iterations large, our methods achieve better losses than the baseline methods.
Loss Iters Dataset: Dialysis
Loss Time
Dataset: Dialysis
10300
10300
10200
10200
losses
10000
10000
9900
9900
0.0
0.5
2.0
2.5
iters
elapsed times (sec)
NewttonExact
NewwtonProx
NewtonQuasi
Quedraticsurrogate (Ourst
Cubiesurrogete (Ours)
Figure 18: Optimization on the Dialysis dataset with A and A2 5.0. Although our methods make less progress toward the minimum loss per iteration (left plot), we are significantly faster than other methods in terms of elapsed time (wall clock) due cheap evaluation cost per iteration. For ease figure reading, we only give partial plot with few iterations. When the number of iterations large, our methods achieve better losses than the baseline methods.
41

---

Loss iters Dataset: Dialysis
vs Time
Dataset: Dialysis
10300
10300
losses
10100
10100
0000000
10300
10000
9900
9900
10
0.0
0.5
1.0
1.5
2.0
25
iters
elapsed times (sec)
NewtonProx
NewtonQuesi
QuatraticSurogstte (Ours)
CubicSurrogate (Ours)
Figure 19: Optimization on the Dialysis dataset with A = 1.0 and X2 = 1.0. The exact Newton method can applied solve the l-rreuularized problems, so only compare with quasi Newton and proximal Newton. Although our methods make less progress toward the minimum loss per iteration (left plot), we are significantly faster than other methods in terms of elapsed time (wall clock) due to cheap evaluation cost per iteration. For ease of figure reading, we only give a partial plot with a few iterations. When the number of iterations is large, our methods achieve better losses than the baseline methods.
Loss VS Iters Dataset: Dialysis
Loss vs Time Dataset: Dialysis
10300
10300
10200
10200
losses
10150
10150
10050
10050
0000
10000
10000
1.0 elapsed times (sec)
2.5
iters
NewtonPrrox
NewtonQuesi
QuadraticSurrgate
CubicSurrogate (Ours)
Figure 20: Optimization on the Dialysis dataset with x = 1.0 and X2 5.0. The exact Newton method can be applied to solve the & l-rrgularized problems, so we only compare with quasi Newton methods make less progress toward the minimum loss per
and proximal Newton. Although iteration (left plot), we are significantly faster than other methods in terms of elapsed time (wall clock) due to cheap evaluation per iteration. For ease of figure reading, only give partial plot with a few iterations. When the number of iterations is large, our methods achieve better losses than the baseline methods.
42

---

D.2
Variable Selection for the CPH Model
D.2.1
Results on Dialysis
Support Size vs train cindex Dataset: Dialysis
Support Size vs test cindex Dataset: Dialysis
0.75
0.75
train cindex
0.70
0.70
test test_cintn cindexiexxx
0000.00
0.65
0.65
0.60
0.60
0.55
Support Size
Support Size
RegimALessoCor
sksurvCoxnet
BzamScorch (Ours)
Figure 21: 5-fold Cross Cross- validation on Dialysis dataset. Comparision with other cox models, metric: CIndex
Support Size train cindex Dataset: Dialysis
Support Size test cindex Dataset: Dialysis
1.0
0.75
0.9
train. cindex
test test_ cindex
0.70
0.8
0.65
0.7
0.60
0.6
10
10
10
10a
10a Support Size
10
Support Size
sisurviree
ASsurVRFF
staurvGaST
sksurVKeiveSVM
sksurvfastSVM
BeamSeerrhh (Ours)
Figure 22: 5-fold cross-validation on Dialysis dataset. Comparision with non-cox models, metric: CIndex
43

---

Support Size train ibs Dataset: Dialysis
Support Size vs test ibs Dataset: Dialysis
0.19
0.19
0.18
0.18
train ibs
test ibs
0.17
10.17
0.16
0.16
0.15
0.15
0.14
0
10
20
30
10
20
30
Support Size
Support Size
segimALessoCoo
aboss
atsurvCoxnet
BzamScarch (Ours)
Figure 23: 5-fold Cross-validation on Dialysis dataset. Comparision with other cox models, metric:
IBS
Support Size vs train ibs Dataset: Dialysis
Support Size vs test ibs Dataset: Dialysis
0.20
0.15
test 139
0.18
train
0.10
0.16
0.05
0.14
101
10
105
10
103
10
Support Size
Support Size
sisurv Tree
sksurvRSF
sksurvGesT
BeamSearch (Durs)
Figure 24: 5-fold Cross-validatio on Dialysis dataset. Comparision with non-cox models, metric: IBS
Support Size vs train_lose Dataset: Dialysis
Support Size vs test loss Dataset: Dialysis
1.65
1.90
1.88
1.60
loss
1.86
test loss
1.55
1.84
1.50
1.82
1.80
1.45
Support Size
Support Size
SogimALossoCox
sksurvCoxnet
BoamSearch
Figure 25: 5-fold Cross-validation on Dialysis dataset. Comparision with other cox models, metric:
CPH Loss
44

---

D.2.2
Results on EmployeeAttrition
Support Size train cindex Dataset: EmployeeAttrition
Support Size test cindexx Dataset: EmployeeAtrittin
0.90
0.90
train, cindex
test_ _s__cindxx
0.85
0.85
0.80
0.80
0.75
10
20
30
10
20
30
Support Size
Support Size
sigimALassoooo
abess
sksurvCoxnet
BeamSearch (Ours)
Figure 26: 5-fold Cross-validation on EmployeeAtrition dataset. Comparision with other cox models,
metric: CIndex
Support Size vs train cindex Dataset: EmployeeAtrritin
train_cindex
Support Size vs test cindex
Dataset: EmployecAtrition
0.95
0.95
train cindex
0.90
test test_cinnee cindex
0.90
0.85
0..55
0.80
0.80
0.75
10
10a
10
10
10a
10
Support Size
Support Size
sksunvWieiveSVM
BearSSerrch (Ours)
Figure 27: 5-fold Cross-validatee EmployeeAntrition dataset. Comparision with non-cox models,
metric: CIndex
45

---

Support Size train ibs Dataset: EmployeeAtrittion
Support Size test ibs Dataset: EmployeeAtritiin
0.200
0.175
0.175
0.150
train ibs
test_ibs
0.150
0.125
0.125
0.100
0.100
0.075
0.075
10
20
30
10
20
30
Support Size
Support Size
segimALessoCoo
aboss
atsurvCoxnet
BzamScarch (Ours)
Figure 28: 5-fold Cross-validation on EmployeeAtrrition dataset. Comparision with other cox models,
metric: IBS
Support Size vs train ibs Dataset: EmployeeAtrition
Support Size test ibs
Dataset: EmployeeAtrrition
0.20
0.15
0.15
test lbs
train
0.10
0.10
0.05
0.05
10
10a
10
10
103 Support Size
10
Support Size
sisurvv Tree
sksurvRSF
sksurvGeST
BeamSearch (Durs)
Figure 29: 5-fold
on EmployeeAttrition dataset. Comparision with non-cox models,
metric: IBS
Support Size vs train loss Dataset: EmployeeAttrition
Support Size vs test loss
Dataset: EmployeeAttrition
1.95
1.90
1.6
train loss
test Joss
1..55
1.5
1.80
1.75
1.4
Support Size
Support Size
SkgimALossoCox
sksurvCoxnet
BoamSearch
Figure 30: 5-fold Cross-validation on EmployeeAtrition dataset. Comparision with other cox models,
metric: CPH Loss
46

---

D.2.3
Results on Kickstarterl
Support Size vs train cindex Dataset: kickstarter1
Support Size vs test cindex Dataset: kickstarter1
1.0
1.0
train cindex
0.9
test_ cindex
0.9
0.8
0.8
0.7
0.7
0.6
0.6
20
80
20
Support Size
Support Size
siglinAeessoCo
abess
sisurvCuxnet
BeanSeerch (Ours)
Figure 31: 5-fold Cross-validation on kickstarterl dataset. Comparision with other cox models,
metric: CIndex
Support Size train cindex Dataset: kickstarter1
Support Size test cindex
Dataset: kickstarterl
1.0
0.95
0.90
train, cindex
test test_ ciindxxxx
0.9
0.85
0.80
0.8
0.75
10a
10a
103
10
10a Support Size
10
Support Size
sksurvTree
SksurvRSF
sksurvGBST
sksurVNeiveSVM
sksurvFastSVM
BearnSarch (Ours)
Figure 32: 5-fold Cross-validation on kickstarterl dataset. Comparision with non-cox models, metric: CIndex
47

---

Support Size train ibs Dataset: kickstarter1
Support Size test ibs Dataset: kickstarter1
0.4
0.4
0.3
0.3
train ibs
test Libbb
0.2
0.2
0.1
0.1
0
20
40
60
80
0
20
40
60
80
Support Size
Support Size
SigimALessoCoo
abass
alsurvCoxnet
BzamScarch (Ours)
Figure 33: 5-fold Cross-validation on kickstarter/ dataset. Comparision with other cox models,
metric: IBS
Support Size vs train ibs Dataset: kickstarter1
Support Size vs test ibs
Dataset: kickstarterl
0.15
0.150
0.125
0.10
test. ibs
0.100
train
0.075
0.05
0.050
0.025
10
103 Support Size
10
10
10a Support Size
10
sisurvires
sksurvRSF
sksurvGesT
BeamSearch (Durs)
Figure 34: 5-fold Cross-validation on kickstarterl dataset. Comparision with non-cox models, metric: IBS
Support Size vs train_loss loss Dataset: kickstarter1
Support Size vs test, loss
Dataset: kickstarterl
4.0
2.8
train loss
3.5
loss
2.6
2.4
3.0
2.2
2.5
2.0
20
80
20
80
Support Size
Support Size
SkgimALossoCox
sksurvCoxnet
BoamSearch
Figure 35: 5-fold Cross-validation on kickstarterl dataset. Comparision with other cox models,
metric: CPH Loss
48

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes],
Justification:
The abstract and intro clearly state the contributions.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Limitations are discussed above the conclusion section.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes].
49

---

Justification:
All proofs are provided in the appendix. All assumptions are provided.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibiliyy
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Experimental setups are clearly stated the Appendix. Code for the method and paper's s experiments is included in public GitHub repository.
Guidelines:
The answer NA means that the paper does include experiments. If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
50

---

Answer: [Yes]
Justification: All baselines and real-world datasets are publicly accessible. The synthetic data generation process is clearly stated in the Appendix. Code for the method and paper's experiments included public GitHub repository.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
Senzidds//ooeSubmiisg
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
7 /mips.cccpublic//u..dde sssp://gguiieeeoeootusuiiiaiiaaa ssinisolicyyyyyc for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/eetails
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
This is described in the experimental section and the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes].
Justification: Error bars are present.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
51

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative emror rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: This reported in the manuscript. See Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics bictps://neurr: http://twrrrpppppppiiietec/////t www.deeeeneeeeeeeeenennnn [Yes]
Justification:
The paper follows the full code of ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction..
10. Broader Impacts
Question: Does the paper discuss both potential positive socictal impacts and negative societal impacts of the work performed?
Answer: [Yes]
The paper involves survival analysis. We have never heard of survival analysis having negative societal impacts. It has positive impacts, particularly in medicine and reliability analysis for equipment/machineery
Justification:
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no socictal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
52

---

The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [Yes]
Justification: The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: The paper develops its own assets. Other works (datasets and baseline software packages) properly cited.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the Datasets.papersstt paperswi thcode. cww/ datasets package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
53

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Question:
Answer:[Yes)
Justification:
Documentation is provided with our code.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/mooell as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA].
Justification:
The paper does not use human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
[NA].
Justification:
This is not relevant to us.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
54


Reviews:
[{"review_id": "2LsdbDfAis", "reviewer": "Reviewer_EDwF", "summary": "The authors propose an optimization of the Cox proportional hazards model based on minimizing surrogate functions obtained from exploiting the Lipschitz continuity property of the first and second order partial derivatives of the loss wrt coefficients. The authors show that the optimization works for sparse penalties or constrained problems like using a cardinality constraint. The authors conduct experiments to show the efficiency and stability of the proposed optimization procedure and empirically show high performance against CoxPH with Newton-Raphson as well as non-linear predictors.", "questions": "- With the proposed method, is it possible to obtain the variance-covariance matrix to construct confidence intervals for the coefficients? If so, the calibration of these would be interesting to look at.\n- Do I understand it correctly that the survival times have to be sorted for the first and second derivative to be O(n) in the current implementation? Would that imply that the complexity of the algorithm is O(n log n)?\n- In future outlook time-varying features are named as an application, does this amount to just changing the risk set in the derivatives as is the case for the loss function?", "limitations": "The authors have addressed some of the limitations of the model.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The size of modern data sets wrt n and p can be a limiting factor for applying the CoxPH model, so innovation in optimization strategy is needed\n- Strongly correlated variables can be an issue for solving the CoxPH model due to colinearity\n- The method works well compared to standard CoxPH empirircally\n- Introducing beamsearch for l0 regularization and covariate selection due to a well behaved loss function leads to interesting performance improvements for large data sets and prevents overfitting", "weaknesses": "- Comparison to an O(n) gradient descent framework like fastCPH (Yang et al) and BigSurvSGD would be helpful in the experiment section (including the lassonet implementation) to substantiate the claim that it is slow (in runtime) compared the proposed procedure and compare stability.\n- When comparing performances wrt survival metrics (c-index, ibs), including the lasso version of the optimization model would be helpful to understand the exact difference between standard lasso cox and the proposed optimization\n- A discussion of the tightness of the bound seems necessary when optimizing the bound instead of the loss function", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for reviewing our paper. Please see below for our answers to your questions.\n\n1. **Comparison with fastCPH and BigSurvSGD on computational efficiency**\n\nThank you for pointing out these two baselines.\nWe will cite them during revision.\n\nRegarding BigSurvSGD, below is a comparison between BigSurvSGD and our methods on Experiment 4.1.\nBigSurvSGD finishes running in a few seconds with low-quality solutions.\nTo compare convergence speed, we set the time limit for our methods so that they run with the same or less amount of time than BigSurvSGD.\nWe can see that our methods achieve much smaller training losses than bigSurvSGD within the same or less wall-clock time, especially using the cubic method.\nWe also want to point out that we implement our methods in python (we loop through all coordinates sequentially) whereas BigSurvSGD is implemented in C++.\n\n| Method      | lambda1 | lambda 2 | train loss | time (s) |\n| ----------- | ----------- | ----------- | ----------- | ----------- |\n| bigSurvSGD | 0 | 1 | 12274 | 18.1 |\n| ours (quadratic) | 0 | 1 | 11054 | 4.2 |\n| ours (quadratic) | 0 | 1 | 10872 | 8.4 |\n| ours (quadratic) | 0 | 1 | 10775 | 16.7 |\n| ours (cubic) | 0 | 1 | 10760 | 4.0 |\n| ours (cubic) | 0 | 1 | 10707 | 7.9 |\n| ours (cubic) | 0 | 1 | 10690 | 15.7 |\n| - | - | - | - | - |\n| bigSurvSGD | 1 | 5 | 10859 | 17.6 |\n| ours (quadratic) | 1 | 5 | 11100 | 4.1 |\n| ours (quadratic) | 1 | 5 | 10935 | 7.6 |\n| ours (quadratic) | 1 | 5 | 10854 | 14.5 |\n| ours (cubic) | 1 | 5 | 10855 | 3.9 |\n| ours (cubic) | 1 | 5 | 10808 | 7.15 |\n| ours (cubic) | 1 | 5 | 10795 | 14.6 |\n\nNote that bigSurvSGD uses stochastic gradient descent.\nStochasticity introduces random noise and does not ensure monotonic loss decrease.\nThis is problematic if we want to use SGD inside beamsearch for the variable selection task in presence of highly correlated features.\nMoreover, for gradient descent, we have mentioned in Line 84-85 that it is hard to pick the right step size, which significantly impacts the convergence speed.\n\nRegarding FastCPH (LassoNetCoxRegressor from the lassonet package), this uses a neural network instead of a pure linear model. \nIt is not allowed to specify no hidden layers for FastCPH (\\# of neural network layers $>$2; see Line 19 in the file model.py in the lassonet package on GitHub).\nTherefore, we cannot use this package to train a linear Cox model.\n    \n2. **Comparison with lasso on solution quality**\n\nIn Section 4.2, we have already included the lasso baseline.\nSee the baseline sksurvCoxnet (the default $\\ell_1$ ratio is $1.0$).\nDuring revision, we will make it clearer that skusrvCoxnet is the lasso baseline.\n    \n3. **Tightness of bounds when optimizing the bound instead of the loss function**\n\nAt each iteration, all current methods (existing Newton methods or our methods) try to minimize a function that approximates the original loss function.\n    \n**Existing Newton methods** They approximate the original loss function through a 2nd-order Taylor expansion (see Line 87-94).\nThe exact Newton method computes the Hessian exactly while quasi and proximal Newton methods compute the Hessian approximately.\n\n**Our quadratic surrogate function** Here, we also make a 2nd-order Taylor expansion but with two differences.\nFirstly, we only do the approximation in a single coordinate.\nSecondly, the expansion's 2nd-order coefficient is an upper bound on the exact 2nd-order derivative.\nFor tightness, if we use the exact 2nd-order derivative instead of the upper bound, we can get a better approximation, which we discuss next.\n\n**Our cubic surrogate function**\nHere, we make a 3rd-order Taylor expansion on a single coordinate.\nWe use the exact 1st- and 2nd -order partial derivatives and also use an upper bound on the 3rd-order partial derivative.\nFor tightness, we obtain a more accurate approximation of the Cox function than our quadratic surrogate function, *at the neighborhood of $x$*.\n\nBoth of our methods ensure monotonic loss decrease and global convergence while existing Newton methods do not.\nPlease see the proof sketch in our general response.\n\n4. **Can we construct confidence intervals for coefficients?**\n\nThis is not the goal of our paper, but yes, it is possible, *after the optimal solution is obtained*.\n\nAfter the optimal solution $\\boldsymbol{\\beta}^*$ is obtained, we first calculate the Hessian $\\boldsymbol{H}$ with respect to the loss function.\nNext, we calculate the standard error vector $\\boldsymbol{s} = \\text{sqrt}\\left[ (\\boldsymbol{H}^{-1}).\\text{diag}() \\right]$, where $\\boldsymbol{Z}.\\text{diag}()$ means taking the diagonal of the square matrix $\\boldsymbol{Z}$ as a vector.\nThe confidence interval is then $\\boldsymbol{\\beta}^* \\pm \\alpha \\boldsymbol{s}$, where $\\alpha \\in \\mathbb{R}_+$ corresponds to different confidence levels.\nPlease see the book In All Likelihood by Yudi Pawitan for details.\nAlso, please see an implementation in the lifelines package on GitHub.\n\n5. **Is sorting needed during optimization?**\n\nPlease look at the our answer in the general response.\n\n6. **How can this work be extended to time-varying features?**\n\nIt takes more than preprocessing and passing a new dataset to our algorithm.\nTime-varying features are not our focus in this work, but we will try our best to answer this question.\nFor details, please see the lifelines documentation page on this topic.\n\nFor time-invariant features, we only care about time duration, and reverse cumulative sum works fine.\nHowever, for time-varying features, we cannot use the vanilla reverse cumulative sum.\nAt different times, new features can enter the loss function (while old features exit).\nOur optimization techniques are still applicable, but we need to write a customized reverse cumulative sum function.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "5Hal4Khqzf", "reviewer": "Reviewer_GRak", "summary": "This paper presents a new optimization algorithm for the Cox model, which is a classical algorithm for survival analysis presented in 1972.", "questions": "N/A", "limitations": "Discussed in page 9.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "This paper is well-written.  I am not an expert on optimization algorithms for convex functions, but I think I could understand the proposed algorithm and I enjoyed reading this paper.\n\nThe proposed algorithm converges faster than other methods due to the new idea of exploiting high order derivatives (up to the third order derivatives) in the Cox model, and it utilizes the surrogate functions for the second and third derivatives.  The proposed algorithm is reasonable and the experimental results are convincing.", "weaknesses": "As a researcher in the field of survival analysis, I would like to note that the impact of this paper is moderate within the field of survival analysis. As mentioned on page 9 of this paper, the Cox model is based on a strong assumption, the proportional hazard assumption, which is unlikely to hold in practice. Therefore, many state-of-the-art models have been developed using much weaker assumptions, and the Cox model is becoming obsolete. However, I would also like to note that the Cox model still shows a certain level of popularity thanks to its simplicity and interpretability.\n\nMisc:\n- Line 52: we are -> ours are\n- Lines 93-94: the formula for quasi Newton method can be simplified by using the diag() function as in the formula for the proximal Newton method.\n- Line 116: result result\n- Line 124: a comma is missing\n- Line 195: wise -> wide\n- Line 257: it is better to clarify which correlation coefficient is used: Pearson, Spearman, or any other one?\n- Page 15: nothing is changed in the statement \"move the partial derivative operator\"\n- Lines 539-540: Redundant vertical lines in the inequality.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for reviewing our paper. Please see below for our answers to your questions.\n\n1. **Utility of Cox Model**\n\n We want to emphasize that we agree with what the reviewer has commented and simply want to continue this conversation. We think there are three aspects that make the Cox model still worthwhile and relevant: \n\n(1) The Cox model is relevant in high stakes domains such as medicine. In these domains, professionals need to interpret the results and understand feature importance. One way to achieve this is to obtain as small number of coefficients as possible without losing predictive performance. In Section 4.2, we have pushed the frontier of this variable selection task.\n\n(2) The Cox model is useful when served as a baseline to compare against more advanced and complicated models. If the more complicated model doesn't outperform the Cox model, we probably want to use the simpler Cox model, according to the principle of Occam's Razor.\n\n(3) Even if we do not deploy the Cox model in the end, it is still beneficial to use during the data exploration stage. When data scientists deal with messy real-world datasets for the first time, they want to perform assessments to see whether there are any abnormalities. It is cumbersome to explore if they use a much more complicated model. Moreover, fixing the abnormalities found by the Cox model could potentially boost the performance of the complicated one.\n \nWe totally respect the reviewer's expertise in survival analysis. We agree with the reviewer that there exist survival models that are not based on the proportional harzard assumption. In Appendeix B Related Work, we have discussed some SOTA survival models. For example, OST and OSST are two decision tree methods. Besides the Cox model, there are also the accelerated failure time (AFT) model and the Aalen’s additive model. Lastly, other models include the ensemble method (SurvivalQuilts) and neural networks (DeepSurv). If you think there are other SOTA survival models we missed, please let us know, we are happy to cite them during the revision.\n   \n2. **Line 93-94, formula for quasi-Newton method**\n\nThank you for this suggestion.\nWe are afraid our current notation cannot be further simplified.\nIn our definition, $\\text{diag}(\\cdot)$ (see Line 94) takes a vector as input and outputs a square matrix whose diagonal equals to this input and other entries should be $0$. However, $\\nabla^2_{\\boldsymbol{\\eta}} \\ell(\\boldsymbol{\\eta})$ is a square matrix.\nThus, we cannot send $\\nabla_{\\boldsymbol{\\eta}}^2 \\ell(\\boldsymbol{\\eta})$ as input to $\\text{diag}(\\cdot)$.\nIn the paper, we had some comments (in green color) to help clarify things.\nThis was the most succinct notation we could come up. If you have another idea to simplify Line 93-94, please let us know, and we are very happy to incorporate that during revision.\n    \n3. **Line 257, definition of feature correlation**\n\nThe correlation used in our work is neither based on Pearson nor Spearman. Please see Line 709-714 in Appendix C2 for the context where the feature correlation is used. For the synthetic dataset, the feature $\\boldsymbol{x}_i$ is sampled from a Gaussian distribution: \n\n$\\boldsymbol{x}_i \\sim \\mathcal{N}(0, \\boldsymbol{\\Sigma})$, \n\nwhere $\\boldsymbol{\\Sigma}$ is the the (Toeplitz) covariance matrix with $\\Sigma_{jl} = \\rho^{|j - l|}$, and $\\rho \\in (0, 1]$ is the correlation parameter.\nWhen $\\rho$ is large, like $\\rho=0.9$ in our paper, the features become highly correlated.\nThis serve as a great benchmark in Section 4.2 to test how well different algorithms can recover the sparse coefficients.\nAs you have seen, our new method significantly outperforms other methods.\n\n4. **Page 15, math typos associated with Appendix A.1.1**\n\nPlease look at our answer in the general response.\n\n5. **Other minor writing typos**\n\n We will fix them during revision. Thank you very much for point these out!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "6i3s70o0WX", "reviewer": "Reviewer_GBdX", "summary": "This paper explores the optimisation of the Cox model. Through careful mathematical analysis, the authors identified efficient ways to calculate the exact derivatives and surrogate loss functions necessary for efficient optimisation, addressing existing strategies' imprecision and time limitations.", "questions": "- It would be beneficial to add a pointer on the proof of convexity (line 82) of the nll.\n- Doesn't the cumulative reverse sum require times to be sorted, i.e. an initial O(n log(n)) computation?\n- Typo Line 195 - 'wise' should be 'wide'.\n\n- In Appendix A1.1.1, the 3rd equation, should have for second element $\\eta_i$ not $\\eta_j$; then no difference between 4th and 5th", "limitations": "As an optimization paper, the authors mention that the work's limitations are the same as those of the standard Cox model.", "rating": 9, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "**Clarity**   \nDespite not being familiar with the literature on optimising the Cox model, the problem is well-described, and the mathematical expressions are rigorous and well-detailed.\n\n**Relevance**\nThe Cox model remains a pillar in multiple fields; ensuring its stability and quick convergence is a significant contribution.", "weaknesses": "Excellent paper, only a few minor typos and questions, see below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for reviewing our paper. Please see below for our answers to your questions.\n\n1. **Proof of convexity of negative log likelihood (nll)** \n\nThe essential part of the nll is the logSumExp function, which is defined as $f(\\boldsymbol{x}) = \\log(\\sum_{i=1}^m \\exp(x_m))$.\nIt is well known the logSumExp function is convex (see Optimization Models and Applications by Laurent El Ghaoui, 2017).\n\nMoreover, the composition of a convex function with an affine (also linear) function is still convex (see the proof on Math Stack Exchange with question ID 654201). Therefore, $\\log(\\sum_{j \\in R_i} e^{\\boldsymbol{x}_j^T \\boldsymbol{\\beta}})$ is convex.\n\nLastly, the sum of two convex functions is still a convex function.\nTherefore, $\\sum_{i=1}^n \\delta_i \\left[ \\log(\\sum_{j \\in R_i} e^{\\boldsymbol{x}_j^T \\boldsymbol{\\beta}}) - \\boldsymbol{x}_i^T \\boldsymbol{\\beta} \\right]$ is convex.\n\n2. **Is sorting needed during optimization?**\n\nPlease look at our answer in the general response.\n\n3. **Page 15, math typos associated with Appendix A.1.1**\n\nPlease look at our answer in the general response.\n\n4. **Minor writing typo**\n\nThank you for your careful reading. We will fix this (wise $\\rightarrow$ wide) during revision.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "v91NlTZeB1", "reviewer": "Reviewer_qYft", "summary": "The authors propose an alternative optimization method for the Cox proportional hazards model. They derive quadratic and cubic upper bounds on the loss and minimize these upper bounds with respect to a single model parameter at a time (similar to coordinate descent) via explicit formulas. They test their method on standard survival analysis benchmarks and also apply it to a feature selection problem with highly correlated features.", "questions": "Is there a mathematical justification for why we should expect the method to be robust to highly correlated features?\n\nPlease also address the concerns listed in the Weaknesses section.", "limitations": "Limitations are briefly discussed (lines 278-282), mostly regarding limitations of the Cox model itself.", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The topic addressed by the paper (survival analysis with large datasets/high dimensionality) is relevant to the community. The prose of the paper was generally easy to follow. The background information provided on survival analysis should also make it easier for a non-expert to understand the paper.", "weaknesses": "I do not believe the paper makes a significant contribution to the literature.\n\nThe runtime reduction results for computing derivatives of the partial likelihood are trivial. The example given on lines 147-148, showing a \"surprising\" improvement from $O(n^2)$ to $O(n)$ computation time, seems to use a deliberately inefficient method for computing derivatives is a point of comparison: there is no need to compute the full Hessian (which is where the authors derive the $O(n^2)$ runtime) when we only want some diagonal terms, and indeed one would expect computing derivatives of a function of the form $\\sum_{i=1}^n f_i(x)$ to take $O(n)$ time provided the derivatives of each function in the summation can be computed in constant time. To that end, the only optimization over a naive implementation is the use of partial sums to avoid repeated calculation of sums with $O(n)$ terms, which is a very standard technique.\n\nThe most important problem is that the experimental results on runtime improvements are misleading. The Flchain dataset is a standard benchmark in survival analysis. It contains 7874 datapoints with 9 features, 4 of which are numeric. After standard 1-hot encoding of the categorical features, the features become 38-dimensional. Requiring 10-20 seconds for convergence (Fig. 1, rightmost plot) on such a modestly-sized dataset does not seem reasonable.\n\nI tested this myself using the scikit-survival (sksurv) library. Sksurv implements Newton's method for fitting the Cox model with an L2 penalty, as can be seen from their source code: https://github.com/sebp/scikit-survival/blob/v0.22.2/sksurv/linear_model/coxph.py#L189\n\nI ran the following code based on the sksurv documentation:\n\n```\nfrom sksurv.datasets import load_flchain\nfrom sksurv.linear_model import CoxPHSurvivalAnalysis\nfrom sksurv.preprocessing import OneHotEncoder\nfrom sklearn.impute import SimpleImputer\n\nX, y = load_flchain() # Load the data\nXt = OneHotEncoder().fit_transform(X) # One-hot encode categorical variables\nimputer = SimpleImputer().fit(Xt) # Dataset contains missing values, impute with mean\nXt = imputer.transform(Xt)\n\ncph = CoxPHSurvivalAnalysis(alpha=1.0, verbose=1) # Fit the model with regularization lambda_2 = 1 as in Fig. 1\ncph.fit(Xt, y)\n```\n\nThe fitting procedure converged for me on my laptop (Apple M2 Pro chip) in a Jupyter notebook in 9 iterations and 0.6s wall-clock time. I also had no problems getting the model fit to converge even with no regularization, contrary to the claim that \"the losses blow up when regularization is weak\" for Newton's method (caption of Fig. 1a). Thus, the apparent gains in performance are most likely due to poor implementations of Newton's method/other baselines.\n\nIn addition to these major issues, there are several other problems:\n1. In equation (6), $\\ell$ had been defined as a function of a $p$-dimensional vector $\\beta$, so it is not clear what is meant by making it a function of the $n$-dimensional vector $\\eta$. The reader can figure out what the authors probably meant, but this should be fixed and defined precisely.\n\n2. The probabilistic interpretations given by the authors for some of the derivative expressions have appeared in very early works on the Cox model, see e.g. Section 2 of [1].\n\n3. Equation (17) is exactly a Newton step, just in a single coordinate and with an upper bound on the second derivative instead of the second derivative itself. It's then not at all clear why we should expect this to be an improvement over Newton's method, especially for the specific cases the authors mention (e.g. high correlations among the features).\n\n4. I'm not sure what is meant by \"the analytical solution to this cubic surrogate function has not been well studied\" (lines 189-190). The cubic formula, which can be used to solve a general cubic equation, has existed for several hundred years (https://en.wikipedia.org/wiki/Cubic_equation).\n\n5. Lines 191-192: The claim that minimizing a convex upper bound of the original function will lead to a decrease in the original function as well is incorrect, barring some important additional assumptions. In particular, I don't believe the authors have justified their claim that the proposed method will \"ensure monotonic loss decrease and global convergence\".\n\n**References:**\n\n[1] Schoenfeld, David. Partial Residuals for the Proportional Hazards Regression Model. Biometrika (1982).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for reviewing our paper. Please see below for our answers to your questions.\n\n1. **Runtime reduction results for computing derivatives are trivial**. \n\nPlease look at Line 146.\nThe result is surprising *especially for computing the exact 2nd order partial derivatives*.\n\nEven for the quasi and proximal Newton methods which make diagonal approximation of Hessian (Line 93-94), the complexity is still high.\nWe **need to calculate the full Hessian** in the space of $\\boldsymbol{\\beta}$ via the formula $\\boldsymbol{X}^T H(\\boldsymbol{\\eta})     \\boldsymbol{X}$, whose complexity is $O(n p^2)$ when $H(\\boldsymbol{\\eta})$ is diagonal.\nPlus, the Newton method requires $O(p^3)$ to solve a linear system via the LDL method. In total, it is $O(np^2 + p^3)$.\n    \nIn contrast, we optimize each coordinate with complexity $O(n)$, and in total $O(np)$ for $p$ coordinates.\nThis reduction is not at all trivial, and we support this claim through empirical results on real-world datasets.\n\n2. **Experimental results on runtime are misleading**\n\n*Your preprocessing of the dataset is different from ours; from running on the wrong preprocessed datasets, one can't conclude our experimental results are wrong*.\nAs written in Line 257-259, we perform binary thresholding to preprocess the data.\nYou didn't do thresholding for continuous variables, which is what makes the variables correlated and challenging.\nIn Appendix C2, our preprocessed Flchain dataset has 333 encoded binary features.\nSince you have only 38 features, you're using a totally different dataset.\n\nThe sksurvCPH you used is also not the exact Newton method.\nThey used a trick on Line 463 of the file coxph.py on the GitHub repo.\nThey take half the exact Newton step size if the loss goes up.\nThis trick does not work in general.\nWhen $\\lambda_2 = 0.001$, sksurvCPH runs over 100 seconds without converging.\nThe loss is 14522.97.\nIn contrast, our method runs under 40 seconds, and the loss is 10659.91 (high precision solution).\n\nWe uploaded to Anonymous GitHub the complete code (existing Newton methods and preprocessing steps) to support our claims.\nDue to the rebuttal policy, we only shared the link with the Area Chair, who can verify results in Section 4.1 and this response by running our codes.\n\n3. **Eq. (6), how do we get from $\\boldsymbol{\\beta}$ to $\\boldsymbol{\\eta}?$**\n\nAs mentioned in Line 89, we use an intermediate variable $\\boldsymbol{\\eta}$ with $\\boldsymbol{\\eta = X \\beta}$.\nBy the chain rule, we have $\\nabla_{\\boldsymbol{\\beta}} \\ell(\\boldsymbol{\\beta}) = \\boldsymbol{X}^T \\nabla_{\\boldsymbol{\\eta}} \\ell(\\boldsymbol{\\eta})$ and $\\nabla_{\\boldsymbol{\\beta}}^2 \\ell(\\boldsymbol{\\beta}) = \\boldsymbol{X}^T \\nabla_{\\boldsymbol{\\eta}}^2 \\ell(\\boldsymbol{\\eta}) \\boldsymbol{X}$.\n\n4. **Probabilistic interpretations are not novel**\n\nThank you for pointing out this paper; it does not overlap with the novel contributions in this submission.\nWe will cite this reference during revision.\nThe reference you cited has only basics.\nThe equation at the bottom of page 1 boils down to the optimality condition when the gradient is equal to 0.\n**Our contributions are substantially different**. We show that\n- The 2nd- and 3rd-order partial derivatives have probabilistic interpretations, which is unexpected.\n- The 2nd- and 3rd-order partial derivatives are in the same formulas as the 2nd- and 3rd-moment calculations.\n- The Lipschitz constants can be computed explicitly for the 1st- and 2nd-order partial derivatives.\n\n5. **Why is Eq.17 an improvement over and different from the existing Newton methods?**\n\nSee the Abstract (Line 7-9), Introduction (Line 33-39), and Preliminaries (Line 95-104). To reiterate:\n\n**Problems with existing Newton methods**\n\na) They all have trouble converging due to the vanishing Hessian.\nThis is a well-known issue for Newton-based methods (see left two plots in Fig 1).\n\nb) Even if they converge, none of them can converge with high precision fast enough (Line 100-104; see right two plots in Fig 1).\n\n**Difference over previous Newton methods**\n\na) We always converge to optimal solutions with high precision. For highly correlated features, other methods could have the loss blow up.\n\nb) Without our formula to calculate $L_2$ in Eq.13, it is impossible to implement Eq.17.\n\nc) Without our formula to calculate $L_3$ in Eq.14, it is impossible to implement for Eq.18.\n\n6. **The cubic equation and formula have existed for several hundred years**\n\nThe cubic equation has almost nothing to do with our cubic surrogate function, despite sharing the same word ``cubic''.\nThe tasks are completely different.\nIn our work, we *minimize* a special polynomial whose highest order is 3.\nOne core step is to *find roots to the second-order polynomial* (see Appendix A4 for details).\nIn contrast, the cubic equation and cubic formula are used to *find the roots to the third-order polynomial*.\nWe never used the cubic formula in our work.\nWe also solve $\\ell_1$-regularized cubic surrogate problem, which has not been studied before.\n    \n7. **Why is the method robust to highly correlated features in variable selection, and is there any math theory?**\n\nAt the core, our methods work because they allow solutions to converge with high precision quickly.\nDuring beam search, having high-precision solutions allows us to compare solutions and pass only the best ones to the next stage of support expansion.\nIf one doesn't get high-precision solutions, some medium- or low-quality solutions can be passed to the next stage, which hurts the overall solution quality.\nFor theories, there are two relevant papers ([13] and [47] in our reference).\nWe can combine their proofs to show our beam search on the Cox model is an approximation algorithm, which means there is some quality guarantee on the final solution.\n\nHopefully we have answered all your questions and cleared all your doubts!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We would like to thank the reviewers for their detailed reviews.
We will use this general response to address some common questions and concerns.

1. **Is sorting needed during optimization?**
Thanks to both Reviewer EDwF and Reviewer GBdX for asking this question.

We never perform sorting during any iterations.
We perform sorting \textbf{only once} at the beginning as a preprocessing step.
When sorting, we can either (1) rearrange the row orders of $\boldsymbol{X}$ and $\boldsymbol{y}$ or (2) record the sorting order and do reverse cumulative sum in this new order.
Thus, the complexity at each iteration is $O(n)$ instead of $O(n \text{log}(n))$. Well-known open-source GitHub packages (scikit-survival, lifelines, skglm, etc) also use sorting as a preprocessing step, so we do not take more preprocessing time than other methods.

2. **Page 15, math typos associated with Appendix A.1.1**
    
Thanks to both Reviewer GRak and Reviewer GBdX for finding the typos.

To fix them, we will delete the 3rd equation with the comment ''move the partial derivative operator'' since the line above this already moves the partial derivative operator inside.
Moreover, from the 3rd to the 8th equations, we will change $-\eta_j \Rightarrow - \eta_i$ and $\frac{\partial}{\partial \eta_{k_1}}(\eta_j) \Rightarrow \frac{\partial}{\partial \eta_{k_1}}(\eta_i)$.
These typos do not affect the end result.
We will correct them during revision.

3. **Monotonic Decrease and Global Convergence**
    
Reviewer qYft asked this question, but the answer here also complements our reply to Reviewer EDwF regarding the tightness of our surrogate functions (bounds).
    
a. **monotonic loss decrease**:
We use the quadratic surrogate function $g_x(\Delta x)$ defined in Equation (15) to illustrate this.
The reasoning also applies to the cubic surrogate function $h_x(\Delta x)$ defined in Equation (16).
From Equation (15), we know two facts:

$$f(x + \Delta x) \leq g_x(\Delta x) \text{ for any } \Delta x, \; \text{ and } f(x) = g_x(0).$$

If we define $\Delta \tilde{x} := \text{argmin}_{\Delta x} g_x(\Delta x)$.
Then we have the inequalities below:

$$f(x + \Delta \tilde{x}) \leq g_x(\Delta \tilde{x}) \leq  g_x(0) = f(x).$$

Thus, we have $f(x + \Delta \tilde{x}) \leq f(x)$.
Since $\Delta \tilde{x}$ is the step we take to minimize $g_x(\cdot)$, our original loss function $f(\cdot)$ will decrease monotonically.
        
b. **global convergence**
We use the Monotone Convergence Theorem (MCT) from basic real analysis to prove convergence.
Let $\{x^t\}$ and $\{f(x^t)\}$ be the sequence of solutions and loss values generated by the iterative procedure, where $t$ is the number of iterations.
From last part, we know that the sequence $\{f(x^t)\}$ is monotonically decreasing.
Furthermore, this sequence is bounded below by $f(x^*)$, where $x^*$ is the optimal solution.
Then, by the MCT, the sequence $\{f(x^t)\}$ converges.

c. **Converging to the optimal value**
We show the sequence $\{f(x^t)\}$ will converge to the optimal value $f(x^*)$.
Because the Cox function $f(\cdot)$ is continuous, convergence of $\{f(x^t)\}$ implies the convergence of $\{x^t\}$.
If we define $\Delta x^t := x^{t+1} - x^{t}$, this implies that  $\Delta x^t \xrightarrow{t \to \infty} 0$.
From Equation (17) or Equation (18), we can  deduce that the first order partial derivatives $f'(x^t)$ converges to 0, for every coordinate.
Note that this is the optimality condition for a convex function.
Thus, we have $\lim_{t \rightarrow \infty} f(x^t) = f(x^*)$.

Meta-review: Cox proportional hazards models are a mainstay of statistical modeling in healthcare and beyond. This paper proposes a simple but important improvement to the stability and speed of optimizing these models, and so represents an important advance in applied data science.
Decision: Accept (poster)