FastSurvival: Hidden Computational Blessings in Training Cox Proportional Hazards Models

Authors: jiachang  liu
rui  zhang
cynthia  rudin
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 17:40:00

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.

Main paper: 
 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 a patient will die [44, 38, 10] . In business, it is 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] . A fundamental tool in analyzing such data is the Cox proportional hazards (CPH) model [8] , a linear model under the assumption that features have a 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, 67, 11, 1] ). Moreover, in contrast to black box models, it is both interpretable and accurate.However, with the advent of larger sample and feature spaces and more complex data, new challenges arise in 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 the CPH model to identify important variables [64, 14] , even in presence of highly correlated features. However, current optimization methods for training the CPH model do not meet these needs. Current algorithms [62, 22, 23, 54] , based on the generic Newton'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 we require the precision 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 a precision-efficiency 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 to 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 (O(n)). Moreover, we show both derivatives are Lipschitz-continuous 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 a quadratic 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 cardinality-constrained CPH problems and produce models that are much sparser than the state-of-the-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 a quadratic and a 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(n)) 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 cardinality-constrained 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 i , t i , δ i } n i=1 , where x ∈ R p is the feature vector with length p, t i ∈ R is the observation time, and δ i ∈ {0, 1} is an 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 as the hazard function in survival analysis. The CPH model predicts the hazard h i (t) for sample i in a semiparametric way [62] . For review of related work, please see Appendix B.EQUATIONwhere h 0 (t) is a baseline hazard function shared by all samples, and β ∈ R p is the parameter of interest. The nice thing about the CPH model is that h 0 (t) cancels out if we look only at the ratio of hazards of sample i vs. all remaining samples at time t i , i.e.,EQUATIONwhere R i := {j | t j ≥ t i } is the set of indices whose observation time is greater than or equal to that of sample i. Such ratios are also called partial likelihoods. To estimate the parameter of interest, β, we maximize the joint partial likelihood of all samples with failure events, which can be written asEQUATIONThis is equivalent to minimizing the negative log partial likelihood [62] , which is defined asEQUATIONThe loss function ℓ(β), 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 R i , for each outer summation index i. Such daunting complexity makes it difficult to employ first-order methods such as gradient descent because we cannot 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 a second-order Taylor expansion:EQUATIONThe function f (∆β) can be minimized by solving a linear system:∆ β = -(∇ 2 β ℓ(β)) -1 ∇ β ℓ(β).To reveal the computational nuances more explicitly, we use an intermediate variable η with η = Xβ. Then we can rewrite the approximation function f (∆β) as:EQUATIONAt each iteration, calculating the Hessian matrix ∇ 2 η ℓ(η) requires O(n 2 ) complexity. Past methods on the CPH model have resorted to various approximation strategies by replacing ∇ 2 η ℓ(η) with H(η) to reduce the computational complexity:1. Exact NewtonH(η) = ∇ 2 η ℓ(η) # no approximation 2. Quasi Newton H(η) ij = [∇ 2 η ℓ(η)] ii if i = j 0 otherwise # ignore off-diagonal terms 3. Proximal Newton H(η) = diag(∇ η ℓ(η) + δ), # diagonal upper bound on ∇ 2 η ℓ(η)where diag(•) constructs a matrix with its diagonal equal to the input vector and other entries equal to 0. 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 method [22] has a local quadratic convergence rate, but each iteration can take a long time. Quasi Newton [62] and proximal Newton [51] 1 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.We want to show that∂ 2 ℓ(β) ∂β 2 l = n i=1 δ i   k∈Ri e η k j∈Ri e ηj X 2 kl - k∈Ri e η k j∈Ri e ηj X kl 2   . ∂ 2 ℓ(β) ∂β 2 l = ∂ 2 ℓ(η) ∂β 2 l = n k2=1 ∂ ∂η k2 ∂ℓ(η) ∂β l ∂η k2 ∂β l # Apply chain rule from calculus = n k2=1 ∂ ∂η k2 n i=1 δ i k∈Ri e η k j∈Ri e ηj X kl - n i=1 δ i X il ∂η k2 ∂β l# plug in ∂ℓ(η) ∂β l from the end result for the first order partial derivative above = n i=1 δ i n k2=1 ∂ ∂η k2 k∈Ri e η k j∈Ri e ηj X kl ∂η k2 ∂β l - n k2=1 ∂ ∂η k2 n i=1 δ i X il ∂η k2 ∂β l # distribute n k2=1 ∂ ∂η k 2 (•)∂η k2 ∂β l # second term is 0 because the input to ∂ ∂η k 2 (•) is a constant = n i=1 δ i n k2=1 k∈Ri ∂ ∂η k2 e η k j∈Ri e ηj X kl ∂η k2 ∂β l # move ∂ ∂η k 2 inside the summation of k = n i=1 δ i n k2=1 k∈Ri j∈Ri e ηj e η k 1 k=k2 -e η k e η k 2 1 y k 2 ≥yi j∈Ri e ηj 2 X kl ∂η k2 ∂β l# evaluate the partial derivative using quotient rule from calculus= n i=1 δ i n k2=1 k∈Ri j∈Ri e ηj e η k 1 k=k2 -e η k e η k 2 1 y k 2 ≥yi j∈Ri e ηj 2 X kl X k2l # ∂η k 2 ∂β l = X k2l = n i=1 δ i n k2=1 k∈Ri j∈Ri e ηj e η k 1 k=k2 -e η k e η k 2 1 y k 2 ≥yi j∈Ri e ηj 2 X kl X k2l # move X k2l inside k∈Ri = n i=1 δ i n k2=1 k∈Ri e η k 1 k=k2 j∈Ri e ηj - e η k e η k 2 1 y k 2 ≥yi j∈Ri e ηj 2 X kl X k2l # divide the fraction into two parts = n i=1 δ i n k2=1 k∈Ri e η k 1 k=k2 j∈Ri e ηj X kl X k2l - e η k e η k 2 1 y k 2 ≥yi j∈Ri e ηj 2 X kl X k2l # distribute X kl X k2l inside (•) = n i=1 δ i k∈Ri n k2=1 e η k 1 k=k2 j∈Ri e ηj X kl X k2l - e η k e η k 2 1 y k 2 ≥yij∈Ri e ηj 2 X kl X k2l # exchange summation orders; sum over k 2 and then k= n i=1 δ i k∈Ri n k2=1 e η k 1 k=k2 j∈Ri e ηj X kl X k2l - n k2=1 e η k e η k 2 1 y k 2 ≥yi j∈Ri e ηj 2 X kl X k2l # distribute n k2=1 inside (•) = n i=1 δ i k∈Ri e η k j∈Ri e ηj X 2 kl - k2∈Ri e η k e η k 2j∈Ri e ηj 2 X kl X k2l # for the first term, only k 2 = k is left because of the indicator 1 k=k2 ; for the second term, the expression can be simplified because R i is a set of indices whose time is greater than or equal toy i = n i=1 δ i k∈Ri e η k j∈Ri e ηj X 2 kl - k∈Ri k2∈Ri e η k e η k 2 j∈Ri e ηj 2 X kl X k2l # distribute k∈Ri inside (•) = n i=1 δ i k∈Ri e η k j∈Ri e ηj X 2 kl - k∈Ri e η k j∈Ri e ηj X kl k2∈Ri e η k 2j∈Ri e ηj X k2l # the double sum of the second term can be rewritten as product of two sums ∂β l inside into the inner summation= n i=1 δ i k∈Ri e η k j∈Ri e ηj X 2 kl - k∈Ri e η k j∈Ri e ηj X kl= n i=1 δ i n k3=1 ∂ ∂η k3 k∈Ri e η k j∈Ri e ηj X 2 kl - k∈Ri e η k j∈Ri e ηj X kl 2 X k3l # ∂η k 3 ∂β l = X k3l = n i=1 δ i n k3=1 ∂ ∂η k3 k∈Ri e η k j∈Ri e ηj X 2 kl - ∂ ∂η k3 k∈Ri e η k j∈Ri e ηj X kl 2 X k3l # distribute ∂ ∂η k 3 inside [•] = n i=1 δ i n k3=1 k∈Ri ∂ ∂η k3 e η k j∈Ri e ηj X 2 kl - ∂ ∂η k3 k∈Ri e η k j∈Ri e ηj X kl 2 X k3l # for the first term, move ∂ ∂η k 3 (•) inside the summation = n i=1 δ i n k3=1 k∈Ri ∂ ∂η k3 e η k j∈Ri e ηj X 2 kl -2 k∈Ri e η k j∈Ri e ηj X kl ∂ ∂η k3 k∈Ri e η kj∈Ri e ηj X kl X k3l # apply product rule from calculus to the second term Proof. Recall that the r-th central moment C r is defined as:= n i=1 δ i n k3=1 k∈Ri ∂ ∂η k3 e η k j∈Ri e ηj X 2 kl -2 k∈Ri e η k j∈Ri e ηj X kl k∈Ri ∂ ∂η k3 e η k j∈Ri e ηj X kl X k3l # for the second term, move ∂ ∂η k 3 (•) inside the summation = n i=1 δ i n k3=1 k∈Ri ∂ ∂η k3 e η k j∈Ri e ηj X 2 kl X k3l -2 k∈Ri e η k j∈Ri e ηj X kl X k3l k∈Ri ∂ ∂η k3 e η k j∈Ri e ηj X kl # distribute X k3l inside {•} = n i=1 δ i n k3=1 k∈Ri e η k j∈Ri e ηj 1 k=k3 - e η k e η k 3 j∈Ri e ηj 2 1 y k 3 ≥yi X 2 kl X k3l -2 k∈Ri e η k j∈Ri e ηj X kl k∈Ri e η k j∈Ri e ηj 1 k=k3 - e η k e η k 3 j∈Ri e ηj 2 1 y k 3 ≥yi X kl X k3l # ∂ ∂η k 3 e η k j∈R i e η j = e η k j∈R i e η j 1 k=k3 -e η k e η k 3 j∈R i e η j 2 1 y k 3 ≥yi = n i=1 δ i n k3=1 k∈Ri e η k j∈Ri e ηj 1 k=k3 - e η k e η k 3 j∈Ri e ηj 2 1 y k 3 ≥yi X 2 kl X k3l -2 k∈Ri e η k j∈Ri e ηj X kl k∈Ri e η k j∈Ri e ηj 1 k=k3 X kl - k∈Ri e η k e η k 3 j∈Ri e ηj 2 1 y k 3 ≥yi X kl X k3l # for the last term, distribute k∈Ri (•)X kl inside [•] = n i=1 δ i n k3=1 k∈Ri e η k j∈Ri e ηj 1 k=k3 X 2 kl X k3l - k∈Ri e η k e η k 3 j∈Ri e ηj 2 1 y k 3 ≥yi X 2 kl X k3l -2 k∈Ri e η k j∈Ri e ηj X kl k∈Ri e η k j∈Ri e ηj 1 k=k3 X kl - k∈Ri e η k e η k 3 j∈Ri e ηj 2 1 y k 3 ≥yi X kl X k3l # for the first term, distribute k∈Ri (•)X 2 kl X k3l inside (•) = n i=1 δ i n k3=1 k∈Ri e η k j∈Ri e ηj 1 k=k3 X 2 kl X k3l - k∈Ri e η k e η k 3 j∈Ri e ηj 2 1 y k 3 ≥yi X 2 kl X k3l -2 k∈Ri e η k j∈Ri e ηj X kl k∈Ri e η k j∈Ri e ηj X kl 1 k=k3 X k3l + 2 k∈Ri e η k j∈Ri e ηj X klC r := k∈Ri e η k j∈Ri e ηj X kl - k1∈Ri e η k 1j∈Ri e ηj X k1l r .We need to show that j∈Ri e ηj X k2l∂ ∂β l C r = C r+1 -r • C 2 • C r-1 . ∂ ∂β l C r = ∂ ∂β l k∈Ri e η k j∈Ri e ηj X kl - k1∈Ri e η k 1 j∈Ri e ηj X k1l r = k∈Ri ∂ ∂β l e η k j∈Ri e ηj X kl - k1∈Ri e η k 1 j∈Ri e ηj X k1l r # move the partial differentail operator ∂ ∂β l (•) inside k∈Ri = k∈Ri ∂ ∂β l e η k j∈Ri e ηj X kl - k1∈Ri e η k 1 j∈Ri e ηj X k1l r + e η k j∈Ri e ηj ∂ ∂β l X kl - k1∈Ri e η k 1The 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 it across all possible time points within the follow-up period of interest. This provides a single score summarizing the model's performance over the entire time range. The lower the IBS score, the better the model.  
new formulas for first, second, and third order partial derivatives: 
 As we have mentioned, the reason for the diagonal approximations of ∇ 2 η ℓ(η) 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 η and focus on the Hessian in the feature space β. 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 a burdensome task, we show that the end result is very elegant and has an intuitive interpretation. Out of complexity comes simplicity. We summarize the relevant results in the first following theorem: Theorem 3.1. For the CPH loss function defined in Equation (4), the first, second, and third order partial derivatives with respect to coordinate l are: 1st order partial derivative:EQUATION2nd order partial derivative:EQUATION)3rd order partial derivative:EQUATIONThe proof can be found in Appendix A. The first, second, and third order partial derivatives all have a probabilistic interpretation. Notice that for any i, the coefficients in front of X kl , X 2 kl , and X 3 kl are nonnegative and sum up to 1, i.e., e η k /( j∈Ri e ηj ) ≥ 0 and k∈Ri [e η k /( j∈Ri e ηj )] = 1. Then, we can regard these coefficients as a discrete probability distribution. Thus, for Equation (8) , the term inside [•] resembles the variance or second order central moment formula:E[X 2 ] -(E[X]) 2 = E[(X -E[X]) 2 ].For Equation ( 9), the term inside [•] resembles the skewness or third order central moment formula:E[X 3 ] + 2(E[X]) 3 -3E[X 2 ]E[X] = E[(X -E[X]) 3 ].One may wonder whether for higher orders (order r ≥ 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 A. Lemma 3.2. Let us define C r to be the r-th central moment withEQUATION)Then we can calculate the partial derivative of C r with respect to β l as:EQUATIONFrom Lemma 3.2, we can see why the connection to central moment does not work for higher order partial derivatives. If r = 2, the second term in Equation ( 11) disappears, i.e., C r-1 = C 1 = 0. Therefore, we get ∂C 2 /∂β l = C 3 . However, for r ≥ 3, C r-1 in general is not zero, so we cannot extrapolate this pattern to higher order partial derivatives.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.From the connections to the second and third central moment, we 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 ∂ℓ(β) ∂β l and ∂ 2 ℓ(β)∂β 2 l are O(n) .This is a surprising result, especially for the second order partial derivatives. The intermediate Hessian, ∇ 2 η ℓ(η), takes O(n 2 ) to compute, so we would expect the second order partial derivative,∂ 2 ℓ(β) ∂β 2 j = e T j X T ∇ 2 η ℓ(η)Xe j , would take O(n 2) to compute as well. Yet, the time complexity is just O(n). We use the first order partial derivative formula, Equation (7) , as an example to explain why this happens. We ignore the second term n i=1 δ i X il because it is just a constant. The first term in Equation ( 7) can be rewritten as:EQUATION)Let us focus on the numerator inside the parenthesis for now. For the entire sequence (i = 1, 2, ..., n) of numerator terms, we can obtain all of them together at the cost of O(n) by performing reverse 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 O(n). Finally, multiplying each ratio with δ i and summing up all these products costs O(n) as well. Therefore, the computational cost to calculate the first order partial derivative is O(n). We can apply the same idea to the second order partial derivative formula, Equation 8, to show that the computational complexity is also O(n). Note that the reverse cumulative summation trick has already been explored in [62] for calculating the diagonal of ∇ 2 η ℓ(η) in the sample space η, but this trick has not been used to calculate the partial derivatives in the feature space β.We will later see how this O(n) time complexity allows us to design a second order optimization method whose evaluation cost per iteration is just as cheap as a first order optimization method. Before we discuss that, let us continue and discuss another computational implication of Theorem 3.1.The connection to the central moment calculation allows us to conclude that the first and second order partial derivatives are Lipschitz-continuous. Moreover, we can calculate these Lipschitz constants explicitly. Recall that for a univariate function f (x), we say that the function is Lipschitzcontinuous [5] if there exists L ≥ 0 such that for any two points in the domain of f , i.e., x, y ∈ D(f ), we have |f (x) -f (y)| ≤ L|x -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, |f ′ (x)| ≤ L for any x ∈ R [5] .Not only can we say that the first and second order partial derivatives are Lipschitz-continuous, 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 (8) , its absolute values are bounded by the following formula:EQUATIONFor the third order partial derivatives in Equation (9) , its absolute values are bounded by the following formula:EQUATION)The availability of these Lipschitz constants suggests that we might construct surrogate functions.We now have all the tools at hand to attack the original optimization problem. For a univariate convex f (x), if we have access to L 2 , the Lipschitz constant for its first order derivative, then we can construct the following quadratic surrogate function g x (•) [53] :EQUATIONIf we have access to L 3 , the Lipschitz constant for its second order derivative, then we can construct the following cubic surrogate function h x (•) [53] :EQUATIONA nice thing about these surrogate functions is that their minimizers can be computed analytically:EQUATIONEQUATIONwhere 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 to 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 f (x) 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 of as coordinate descent-type methods. We anticipate these core ideas can be applied to solve a wide range of problems related to the CPH model. In the next subsection, we showcase two problems our algorithms can tackle.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 [15] , MCP [68] , etc. For the ℓ 1 -regularized problems, we can in fact find analytical solutions2 .For the quadratic surrogate function, solving the ℓ 1 -regularized problem in Equation ( 15) is equivalent to solving the following optimization problem (with a = f ′ (x), b = L 2 , and c = x),EQUATION)The solution for the above problem isEQUATIONFor the cubic surrogate function, solving the ℓ 1 -regularized problem of Equation ( 16) is equivalent to solving the following optimization problem (with a = f ′ (x), b = f ′′ (x), c = L 3 , and d = x):EQUATIONwhose solution is:EQUATIONEquation ( 20) is well known in a slightly different format. Equation (22) has not been well studied in the past. We provide derivations for both in Appendix A.The second problem is the cardinality-constrained CPH problem. Recently, the beam search framework (a combination of the beam search method [66] from natural language processing and generalized orthogonal matching pursuit [13] ) has shown promise in finding near optimal solutions for a class of l 0 -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 cardinality-constrained problems, it cannot be applied directly to the CPH model without our coordinate descent methods to select features, especially in the highly correlated settings. 
experiments: 
 We test the effectiveness of our optimization methods on both synthetic and real-world datasets. We run experiments for both regularized and constrained problems mentioned in Section 3.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-the-arts methods, especially for challenging scenarios where features are highly correlated?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 2: exact Newton method, the quasi Newton method, and the proximal Newton method. We run on both ℓ 2 -regularized CPH problems and ℓ 1 + ℓ 2 -regularized CPH problems. The choices of these regularizations are:λ 2 = {0, 1} and λ 2 = {1, 5}.The coefficients are all initialized to be 0. In the main paper, we show results on the Flchain dataset in Figure 1 . More results on other datasets can be found in Appendix D.During each iteration, the baseline methods [62, 51] optimize all coefficients at once, whereas our 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-type baselines sometimes have losses that blow up or increase during the initial phase of optimization. This is a common problem of Newton'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 vs. overall time elapsed, we can see that our 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 O(n)), as we have explained in Section 3.2.We compare our method with both Cox-based methods and other model classes. For Cox-based 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 in all datasets. Synthetic datasets are generated with high correlation level, ρ = 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 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 on other datasets.Figure 2 : Variable selection on synthetic datasets with high correlation (correlation level ρ = 0.9).From left to right, the sample sizes are 1200, 1000, and 800, respectively. The F1 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 F1 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 F1 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 baselines. 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 we 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 accuracy-sparsity tradeoff. We are able to obtain solutions with the smallest number of coefficients without losing predictive performance.Figure 3 : Variable selection on the Employee 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.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 is 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 to train the Cox proportional hazards (CPH) model by constructing 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, it will be interesting 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-continuous 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 [16] , stratifications [40] , and feature interactions [45] . 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: This was reported in the manuscript. See Appendix C.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 runs 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 make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? [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 circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes]Justification: 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/machinery.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or 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.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a 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 is 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 to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: 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)?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 a 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 a 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) are 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 or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.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 package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a 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. 
 
 Appendix: A.1 First, Second, and Third Order Partial DerivativesThe 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. The derivation for the third order partial derivative is shown in Section A.1.3.We want to show that∂ℓ(β) ∂β l = n i=1 δ i k∈Ri e η k j∈Ri e ηj X kl - n i=1 δ i X il . ∂ℓ(β) ∂β l = ∂ℓ(η) ∂β l = n k1=1 ∂ℓ(η) ∂η k1 ∂η k1 ∂β l # apply chain rule from calculus = n k1=1 ∂ ∂η k1 n i=1 δ i log( j∈Ri e ηj ) -η i ∂η k1 ∂β l # plug in the formula for ℓ(η) = n i=1 δ i n k1=1 ∂ ∂η k1 log( j∈Ri e ηj ) -η i ∂η k1 ∂β l# exchange the summation orders; sum over k 1 first and sum over i later= n i=1 δ i n k1=1 ∂ ∂η k1 log( j∈Ri e ηj ) -η i X k1l# evaluate the partial derivative∂η k 1 ∂β l = X k1l = n i=1 δ i n k1=1 ∂ ∂η k1 (log( j∈Ri e ηj )) - ∂ ∂η k1 (η i ) X k1l # distribute ∂ ∂η k 1 (•) inside [•] = n i=1 δ i n k1=1 1 j∈Ri e ηj e η k 1 1 y k 1 ≥yi -1 k1=i X k1l # ∂ ∂η k 1 (log( j∈Ri e ηj )) = 1 j∈R i e η j e η k 1 1 y k 1 ≥yi and ∂ ∂η k 1 (η i ) = 1 k1=i = n i=1 δ i n k1=1 1 j∈Ri e ηj e η k 1 1 y k 1 ≥yi X k1l -1 k1=i X k1l # distribute X k1l inside [•] = n i=1 δ i n k1=1 1 j∈Ri e ηj e η k 1 1 y k 1 ≥yi X k1l - n k1=1 1 k1=i X k1l # distribute n k1=1 inside (•) = n i=1 δ i n k1=1 1 y k 1 ≥yi 1 j∈Ri e ηj e η k 1 X k1l - n k1=1 1 k1=i X k1l # move 1 y k 1 ≥yi to the front = n i=1 δ i k1∈Ri 1 j∈Ri e ηj e η k 1 X k1l -X il# evaluate the two summations over k 1 , i.e., n k1=1 ; the first summation can be simplified because R i is a set of indices whose time is greater than or equal toy i = n i=1 δ i k1∈Ri 1 j∈Ri e ηj e η k 1 X k1l - n i=1 δ i X il # distribute n i=1 δ i inside (•) = n i=1 δ i k1∈Ri e η k 1 j∈Ri e ηj X k1l - n i=1 δ i X il # move e η k 1 to the numerator = n i=1 δ i k∈Ri e η k j∈Ri e ηj X kl - n i=1 δ i X il # change notation k 1 to k A.1.Let us now plug this result back into the original problem:∂ ∂β l (C r ) = k∈Ri ∂ ∂β l e η k j∈Ri e ηj X kl - k1∈Ri e η k 1 j∈Ri e ηj X k1l r -r k1∈Ri ∂ ∂β l e η k 1j∈Ri e ηj X k1l C r-1 # pick up from where we left for the original problem= k∈Ri e η k j∈Ri e ηj X kl - k2∈Ri e η k 2 j∈Ri e ηj X k2l X kl - k1∈Ri e η k 1 j∈Ri e ηj X k1l r -r k1∈Ri e η k 1 j∈Ri e ηj X k1l - k2∈Ri e η k 2j∈Ri e ηj X k2l X k1l C r-1 # plug in the solution to the subproblem= k∈Ri   e η k j∈Ri e ηj X kl - k1∈Ri e η k 1 j∈Ri e ηj X k1l r+1   -r k1∈Ri e η k 1 j∈Ri e ηj X k1l - k2∈Ri e η k 2 j∈Ri e ηj X k2l X k1l C r-1# for the first term, change k 2 into k 1 because both are dummy variables that are independent from each other= k∈Ri   e η k j∈Ri e ηj X kl - k1∈Ri e η k 1 j∈Ri e ηj X k1l r+1   -r k1∈Ri e η k 1 j∈Ri e ηj X 2 k1l - k2∈Ri e η k 2 j∈Ri e ηj X k1l X k2l C r-1# for the second term, moveX k1l inside (•) = k∈Ri   e η k j∈Ri e ηj X kl - k1∈Ri e η k 1 j∈Ri e ηj X k1l r+1   -r    k1∈Ri e η k 1 j∈Ri e ηj X 2 k1l - k1∈Ri k2∈Ri e η k 1 e η k 2 j∈Ri e ηj 2 X k1l X k2l    C r-1 # for the second term, distribute k1∈Ri e η k 1 j∈R i into the two terms inside (•) = k∈Ri   e η k j∈Ri e ηj X kl - k1∈Ri e η k 1 j∈Ri e ηj X k1l r+1   -r   k1∈Ri e η k 1 j∈Ri e ηj X 2 k1l - k1∈Ri e η k 1 j∈Ri e ηj X k1l 2   C r-1# in the second term, k 1 and k 2 are independent dummy variables, so we can turn a double sum of products into a products of sums; we can further simplify this into a square of a sum because the two terms equal to the same value= C r+1 -r • C 2 • C r-1# simplify by using the definition of central momentProof. To show that the first and second order partial derivatives are Lipschitz-continuous, we need to show that the second and third order partial derivatives are bounded, respectively.First order partial derivative is Lipschitz-continuous 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 to each coordinate is bounded. Recall that the second order partial derivative with respect to the l-th coordinate can be expressed as:EQUATIONIt suffices to show that each term inside the bracket is bounded.If we interpret the expression probabilistically, then the coefficientse η kj∈R i e η j in front of X 2 kl and X kl can be thought of as the probability of 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 a k = e η k j∈R i e η j . Then we can rewrite each term inside [•] as:k∈Ri e η k j∈Ri e ηj X 2 kl - k∈Ri e η k j∈Ri e ηj X kl 2 = k∈Ri a k X 2 kl -( k∈Ri a k X kl ) 2The right-hand side is nothing but the variance of {X kl } k∈Ri with respect to the distribution {a k } k∈Ri . Since the variance is always non-negative, we havek∈Ri a k X 2 kl -( k∈Ri a k X kl ) 2 ≥ 0Let us now denote a := min k∈Ri X kl as the minimum of this given set, b := max k∈Ri X kl as the maximum of this given set, Z as a random variable with values restricted to [a, b]. We are going to show thatEQUATION)We achieve this through two steps.First EQUATIONwhere 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 [•], so taking max(•) produces the inequality above.Next, if we drop the assumption that the random variable Z can only take |R i | number of values, then the maximum variance we can achieve is no smaller than before. Mathematically, this means that maxEQUATION)Combining Inequality (24) and Inequality (25) , we arrive at Inequality (23) . Lastly, note that the Popoviciu Inequality [58] tells us that for a bounded random variable restricted to [a, b], the maximum variance it can achieve is (a-b) 2, i.e.,max Z E[Z 2 ] -(E[Z]) 2 ≤ (b -a) 2 4 .This allows us to conclude thatk∈Ri a k X 2 kl - k∈Ri a k X kl 2 ≤ (b -a) 2 4 = (max k∈Ri X kl -min k∈Ri X kl ) 2 4 .Plugging in this inequality into the expression of the second order partial derivative, we can show that the second order partial derivative is bounded:0 ≤ ∂ 2 ℓ(β) ∂β 2 l = n i=1 δ i k1∈Ri e η k 1 j∈Ri e ηj X 2 k1l - k1∈Ri e η k 1 j∈Ri e ηj X k1l 2 ≤ 1 4 n i=1 δ i (max k∈Ri X kl -min k∈Ri X kl ) 2Second order partial derivative is Lipschitz-continuous This proof is similar to the first part above. We need to show that the third order partial derivative with respect to each coordinate is bounded. Recall that the third order partial derivative can be expressed as: ∂ 3 ℓ(β) ∂β 3 l = n i=1 δ i k∈Ri e η ke η k j∈Ri e ηj X kl = k∈Ri a k X 3 kl + 2( k∈Ri a k X kl ) 3 -3( k∈Ri a k X 2 kl )( k∈Ri a k X kl )Similarly, as in the first part, let us now denote a := min k∈Ri X kl as the minimum of this given set, b := max k∈Ri X kl as the maximum of this given set, Z as a random variable with values restricted to [a, b]. Using the exact same logic, we have| k∈Ri a k X 3 kl + 2( k∈Ri a k X kl ) 3 -3( k∈Ri a k X 2 kl )( k∈Ri a k X kl )| ≤ max Z |E[Z 3 ] + 2E[Z] 3 -3E[Z 2 ]E[Z]||| = max Z |E[(Z -E[Z]) 3 ]|The expression inside |•| on 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 a bounded variable.According to [61] , we have the following inequality involving the second and third central moment:E[(Z -E[Z]) 2 ] + E[(Z -E[Z]) 3 ] 2E[(Z -E[Z]) 2 ] 2 ≤ 1 4 (b -a) 2From this, we can derive an upper bound on the third central moment:|E[(Z -E[Z]) 3 ]| ≤ 2E[(Z -E[Z]) 2 ] 1 4 (b -a) 2 -E[(Z -E[Z]) 2 ]For notational convenience, let us denoteV := E[(Z -E[Z]) 2 ], then the right-hand side above can be expressed as a function of V:f (V ) := 2V 1 4 (b -a) 2 -V = (b -a) 2 V 2 -4V 3 Because V is the variance, V ∈ [0, 1 4 (b -a) 2 ]. Additionally, because• is monotonically increasing, the maximum is achieved either at the points where the first order derivative of (b -a) 2 V 2 -4V 3 with respect to V is 0 or at the boundary, 0 and 1 4 (b -a) 2 . Let us calculate the points where the first order derivative is 0:d dV (b -a) 2 V 2 -4V 3 = 2(b -a) 2 V -12V 2 = 0 ⇐⇒ V = (b -a) 2 6 or V = 0To obtain the maximum value achievable, we calculate the values at points V = 0, V = (b-a) 2, and V = 1 4 (b -a) 2 and pick the maximum value afterwards:f (0) = 2 × 0 1 4 (b -a) 2 -0 = 0 f ( 1 6 (b -a) 2 ) = 2 × 1 6 (b -a) 2 1 4 (b -a) 2 - 1 6 (b -a) 2 = 1 6 √ 3 |b -a| 3 f ( 1 4 (b -a) 2 ) = 2 × 1 4 (b -a) 2 1 4 (b -a) 2 - 1 4 (b -a) 2 = 0Therefore, the maximum value achievable is 1 6 √ 3 |b -a| 3 . We now show that upper bound on the absolute value of the third central moment is actually tight by providing with a concrete example. For a random variable Z, letP[Z = a] = 1 4 , P[Z = b] = 1 4, andP[Z = a+b 2 ] = 1 2 . We can verify that E[(Z -E[Z]) 3 ] = 1 6 √3 |b -a| 3 , thus proving that this upper bound is indeed tight. This helps us to arrive at the following inequality:| k∈Ri a k X 3 kl + 2( k∈Ri a k X kl ) 3 -3( k∈Ri a k X 2 kl )( k∈Ri a k X kl )| ≤ 1 6 √ 3 |max k∈Ri X kl -min k∈Ri X kl | 3Therefore, for the upper bound of the third order partial derivative, we have the following explicit formula:| ∂ 3 ℓ(β) ∂β 3 l | ≤ 1 6 √ 3 n i=1 δ i |max k∈Ri X kl -min k∈Ri X kl | 3Let f (x) be a convex function whose first, second, and third derivatives all exist. Let h x (∆x) := f (x) + f ′ (x)∆x + 1 2 f ′′ (x)∆x 2 + 1 6 L 3 |∆x| 3 be the cubic surrogate function [53] of f (x), where L 3 > 0 is the Lipschitz-constant of the second derivative f ′′ (x). Then, the minimum of this surrogate function is achieved at the following point:EQUATIONProof. We discuss three cases: f ′ (x) > 0, f ′ (x) < 0, and f ′ (x) = 0.Case 1 f ′ (x) > 0.If f ′ (x) > 0, then ∆x < 0. For the sake of contradiction, suppose ∆x > 0, which means h x (∆x) achieves its minimum at some point with ∆x > 0. However, we arrive at a contradiction becauseh x (0) = f (x) < f (x) + f ′ (x)∆x + 1 2 f ′′ (x)∆x 2 + 1 6 L 3 |∆x| 3 = h x (∆x), for ∆x > 0.Therefore, the minimum is achieved either at ∆x = 0 or ∆x < 0. However, since ∆x 2 and |∆x| 3 grow slower than |∆x| when ∆x is close to 0, there exists some ∆x < 0 such that f ′ (x)∆x + 1 2 f ′′ (x)∆x 2 + 1 6 L 3 |∆x| 3 < 0. Thus, h x (0) cannot be the minimum value, and we are left with the minimum value achieved at some ∆x < 0. 3 . Note that the second order derivative of h x (∆x) is greater than or equal to 0 sinceIf ∆x < 0, h x (∆x) = f (x) + f ′ (x)∆x + 1 2 f ′′ (x)∆x 2 - 1 6 L 3 (∆x)d 2 d∆x 2 h x (∆x) = f ′′ (x) -L 3 ∆x ≥ 0. # f ′′ (x) ≥ 0 because f (x) is convex, 0 ≤ |f ′′′ (x)| ≤ L 3 ,Therefore, h x (∆x) is a convex function with respect to ∆x when ∆x < 0, and its minimum value is achieved when the first order derivative is 0. When the derivative with respect to ∆x is 0, we havef ′ (x) + f ′′ (x)∆x - 1 2 (∆x) 2 = 0 ⇐⇒ ∆x = f ′′ (x) ± (f ′′ (x)) 2 + 2L 3 f ′ (x) L 3 .Since f ′ (x) > 0, only one rootf ′′ (x)- √ (f ′′ (x)) 2 +2L3f ′ (x) L3satisfing the condition ∆x < 0. Thus, when f ′ (x) < 0, we have∆x = f ′′ (x) -(f ′′ (x)) 2 + 2L 3 f ′ (x) L 3 Case 2 f ′ (x) < 0.If f ′ (x) < 0, we have ∆x > 0 using the same logic as above. 3 . Note that the second order derivative of h x (∆x) is also greater than or equal to 0 sinceIf ∆x > 0, h x (∆x) = f (x) + f ′ (x)∆x + 1 2 f ′′ (x)∆x 2 + 1 6 L 3 (∆x)d 2 d∆x 2 h x (∆x) = f ′′ (x) + L 3 ∆x ≥ 0. # f ′′ (x) ≥ 0 because f (x) is convex, L 3 > 0, and ∆x > 0Therefore, h x (∆x) is a convex function with respect to ∆x when ∆x > 0, and its minimum value is achieved when the first derivative is 0. When the derivative with respect to ∆x is 0, we havef ′ (x) + f ′′ (x)∆x + 1 2 (∆x) 2 = 0 ⇐⇒ ∆x = -f ′′ (x) ± (f ′′ (x)) 2 -2L 3 f ′ (x) L 3 .Since f ′ (x) < 0, only one root-f ′′ (x)+ √ (f ′′ (x)) 2 -2L3f ′ (x) L3satisfing the condition ∆x > 0. Thus, when f ′ (x) < 0, we have∆x = -f ′′ (x) + (f ′′ (x)) 2 -2L 3 f ′ (x) L 3 Case 3 f ′ (x) = 0.When f ′ (x) = 0, the minimum value of h x (∆x) is achieved at ∆x = 0.The explicit formulas for the three cases above can be unified into one succinct formula below:∆x = sgn(f ′ (x)) • f ′′ (x) -(f ′′ (x)) 2 + 2L 3 |f ′ (x)| L 3A.5 Analytical Solution to the ℓ 1 -regularized Quadratic and Cubic Surrogate Problems ℓ 1 -regularized quadratic surrogate problem We have the following ℓ 1 -regularized quadratic surrogate problem:∆x = argmin ∆x a∆x + 1 2 b∆x 2 + λ 1 |c + ∆x|.The solution for the above problem is∆x =    -(a -λ 1 )/b if bc -a < -λ 1 -(a + λ 1 )/b if bc -a > λ 1 -c otherwise .Proof. Since the function a∆x + 1 2 b∆x 2 + λ 1 |c + ∆x| is convex, the condition for this function to achieve the minimum value is for its differential to include 0. The differential of this function is:∂ ∆x a∆x + 1 2 b∆x 2 + λ 1 |c + ∆x| =    a + b∆x + λ 1 if ∆x > -c a + b∆x -λ 1 if ∆x < -c a + b∆x + [-λ 1 , λ 1 ] if ∆x = -c 1.For the first condition, if the differential contains 0, we havea + b∆x + λ 1 = 0 ⇒ ∆x = - a + λ 1 bHowever, because we require ∆x > -c, we have- a + λ 1 b > -c ⇒ bc -a > λ 12. For the second condition, if the differential contains 0, we havea + b∆x -λ 1 = 0 ⇒ ∆x = - a -λ 1 b However, because we require ∆x < -c, we have - a -λ 1 b < -c ⇒ bc -a < -λ 13. For the third condition, if the differential contains 0, we have0 ∈ a + b∆x + [-λ 1 , λ 1 ] ⇒ a + b∆x -λ 1 ≤ 0 ≤ a + b∆x + λ 1However, because we require ∆x = -c, we have -λ 1 ≤ bc -a ≤ λ 1 ℓ 1 -regularized cubic surrogate problem We have the following ℓ 1 -regularized cubic surrogate problem:EQUATIONThe solution to the above problem is∆x =              sgn(d) -b + b 2 -2c(sgn(d)a + λ 1 ) /c if sgn(d)a + λ 1 ≤ 0 sgn(d) b + b 2 + 2c(sgn(d)a -λ 1 ) /c if sgn(d)(a -bd) -1 2 cd 2 > λ 1 sgn(d) b + b 2 + 2c(sgn(d)a + λ 1 ) /c if sgn(d)(a -bd) -1 2 cd 2 < -λ 1 -d otherwise .Proof. Like the first part, since the function a∆x + 1 2 b∆x 2 + 1 6 c|x| 3 + λ 1 |d + ∆x| is convex, the condition for this function to achieve the minimum value is for its differential to include 0. We discuss the differential of this function in two cases: d ≥ 0 and d < 0.• When d ≥ 0, the differential of this function is:∂ ∆x a∆x + 1 2 b∆x 2 + 1 6 c|x| 3 + λ 1 |d + ∆x| =              a + b∆x + 1 2 c∆x 2 + λ 1 if ∆x > 0 a + b∆x + 1 2 c[-∆x 2 , ∆x 2 ] + λ 1 if ∆x = 0 a + b∆x -1 2 c∆x 2 + λ 1 if -d < ∆x < 0 a + b∆x -1 2 c∆x 2 + [-λ 1 , λ 1 ] if ∆x = -d a + b∆x -1 2 c∆x 2 -λ 1 if ∆x < -dWe discuss these 5 cases one by one.1. For the first condition, if the differential contains 0, we havea + b∆x + 1 2 c∆x 2 + λ 1 = 0 ⇒ ∆x = -b ± b 2 -2c(a + λ 1 )c However, because we require ∆x > 0, we haveλ 1 = -(a + b∆x + 1 2 c∆x 2 ) < -aThis means that we can only have one root because the other root violates ∆x > 0:∆x = -b + b 2 -2c(a + λ 1 ) c 2.For the second condition, if the differential contains 0, we have0 ∈ a + b∆x + 1 2 c[-∆x 2 , ∆x 2 ] + λ 1 ⇒ a + b∆x - 1 2 c∆x 2 + λ 1 ≤ 0 ≤ a + b∆x + 1 2 c∆x 2 + λ 1However, because we require ∆x = 0, we have a + λ 1 = 0 3. For the third condition, if the differential contains 0, we havea + b∆x - 1 2 c∆x 2 + λ 1 = 0 ⇒ ∆x = -b ± b 2 + 2c(a + λ 1 )-c However, because we require -d < ∆x < 0, we haveλ 1 = -a + b∆x - 1 2 c∆x 2 > -a ⇒ a + λ 1 > 0This means that we can only have one root because the other root violates the condition-d < ∆x < 0: ∆x = b -b 2 + 2c(a + λ 1 )c Moreover, since the root is between -d and 0, and the coefficients in front of ∆x, -1 2 c, is negative, we have a -bd -1 2 cd 2 + λ 1 < 0. 4. For the fourth condition, if the differential contains 0, we have0 ∈ a + b∆x - 1 2 c∆x 2 + [-λ 1 , λ 1 ] ⇒ a + b∆x - 1 2 c∆x 2 -λ 1 ≤ 0 ≤ a + b∆x - 1 2 c∆x 2 + λ 1However, because we require ∆x = -d, we havea -bd - 1 2 cd 2 -λ 1 ≤ 0 ≤ a -bd - 1 2 cd 2 -λ 1 ⇒ |a -bd - 1 2 cd 2 | ≤ 0 5.For the fifth condition, if the differential contains 0, we havea + b∆x - 1 2 c∆x 2 -λ 1 = 0 ⇒ ∆x = -b ± b 2 + 2c(a -λ 1 )-c However, because we require ∆x < -d, we haveλ 1 = a + b∆x - 1 2 c∆x 2 < a ⇒ a -λ 1 > 0This means that we can only have one root because the other root violates the condition ∆x < -d ≤ 0:∆x = b -b 2 + 2c(a -λ 1 )c Moreover, because the left root is less than -d and the coefficient in front ∆x 2 , -1 2 c, is negative, we have a -bd -1 2 cd 2 -λ 1 > 0.• When d < 0, the differential of this function is∂ ∆x a∆x + 1 2 b∆x 2 + 1 6 c|x| 3 + λ 1 |d + ∆x| =              a + b∆x + 1 2 c∆x 2 + λ 1 if ∆x > -d a + b∆x + 1 2 c∆x 2 + [-λ 1 , λ 1 ] if ∆x = -d a + b∆x + 1 2 c∆x 2 -λ 1 if 0 < ∆x < -d a + b∆x + 1 2 c[-∆x 2 , ∆x 2 ] -λ 1 if ∆x = 0 a + b∆x -1 2 c∆x 2 -λ 1 if ∆x < 0Similar to the previous part when d ≥ 0, we discuss the 5 cases one by one but omit the details because the logic and the reasoning process are exactly the same:1. For the first condition, if the differential contains 0, we have∆x = -b + b 2 -2c(a + λ 1 ) c and a -bd + 1 2 cd 2 + λ 1 < 0 2.For the second condition, if the differential contains 0, we have∆x = -d and |a -bd + 1 2 cd 2 | ≤ λ 13. For the third condition, if the differential contains 0, we have∆x = -b + b 2 -2c(a -λ 1 ) c 4.For the fourth condition, if the differential contains 0, we have ∆x = 0 and a -λ 1 = 0 5. For the fifth condition, if the differential contains 0, we have∆x = b -b 2 + 2c(a -λ 1 ) c and a -λ 1 > 0We can combine the two situations where d ≥ 0 and d < 0 and obtain a unified formula:∆x =                  sgn(d) -b + b 2 -2c(sgn(d)a + λ 1 ) /c if sgn(d)a + λ 1 < 0 0 if sgn(d)a + λ 1 = 0 sgn(d) b + b 2 + 2c(sgn(d)a -λ 1 ) /c if sgn(d)(a -bd) -1 2 cd 2 > λ 1 sgn(d) b + b 2 + 2c(sgn(d)a + λ 1 ) /c if sgn(d)(a -bd) -1 2 cd 2 < -λ 1 -d otherwise .The first and second equations above can be further unified into just one equation sgn(d) -b + b 2 -2c(sgn(d)a + λ 1 ) /c if sgn(d)a + λ 1 ≤ 0, so we finally have ∆x =              sgn(d) -b + b 2 -2c(sgn(d)a + λ 1 ) /c if sgn(d)a + λ 1 ≤ 0 sgn(d) b + b 2 + 2c(sgn(d)a -λ 1 ) /c if sgn(d)(a -bd) -1 2 cd 2 > λ 1 sgn(d) b + b 2 + 2c(sgn(d)a + λ 1 ) /c if sgn(d)(a -bd) -1 2 cd 2 < -λ 1 -d otherwise .as possible without losing accuracy. In this work and especially the variable selection experiments, we push the limit of sparsity-accuracy tradeoff curve for this model class.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.We have a summary of datasets for experiments in Table 1 . The synthetic data used in the paper is generated according to the following process, similar to [71] :Firstly, from a Gaussian distribution N (0, Σ) where the first entry is the mean and the second entry is the covariance matrix with size p × p, we sample features:EQUATIONThe covariance matrix is defined as Σ jl = ρ |j-l| , where ρ ∈ (0, 1] is the correlation parameter. When ρ is large, the features in x i are more correlated. We create a k-sparse coefficient vector β * ∈ R p . The entries of β * are either 1 or 0. If jmod(p/k) = 0, then β * j = 1; otherwise, β * j = 0. Secondly, we generate the death time t i according to the following equation:EQUATIONwhere V i ∼ U (0, 1) (samples are drawn from a uniform distribution on the interval [0, 1]) and s is a hyperparameter. In 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 a uniform distribution: C i ∼ U (0, 1). If the death time is bigger than the censoring time, we have the indicator equal to 1; otherwise, we have the indicator equal to 0. Specifically, we have:EQUATIONAfterwards, we change the death time to observation time, taking into consideration of censoring:EQUATIONWe form a triplet (x i , t i , δ i ) and return this triplet as one sample.Real-world survival data:• Flchain: Use of nonclonal serum immunoglobulin free light chains to predict overall survival in the general population [10] . The event is death. • Kickstarter1: Data from a popular crowdfunding platform, used to predict project success [46] . We used the version from https://dmkd.cs.vt.edu/projects/survival/ data/.Collection and Setup: We ran 5-fold cross-validation (random seed 0) on the following datasets: Dialysis, Flchain, Kickstarter1, EmployeeAttrition, SyntheticHighCorrHighDim1, SyntheticHighCor-rHighDim2, SyntheticHighCorrHighDim3. In order to create highly correlated features, we encoded 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:• Abess: 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 primary_model_f it_max_iter to be 20, approximate_N ewton to 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 l1_ratio to be 1.0, alpha_min_ratio to be 0.01. All other parameters were set to the default. • SkglmALassoCox: We ran this algorithm with 9 different L1 regularization penalty parameters (alpha): 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100. All other parameters were set to the default. • SksurvTree: We ran this algorithm with 8 different configurations: max depth limit, d, ranging from 2 to 9, and a corresponding maximum leaf limit 2 d . The random state was set to 2024 and all other parameters were set to the default. • SksurvRSF, SksurvGBST: We ran this algorithm with 8 × 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. • SksurvNaiveSVM, SksurvFastSVM: We ran this algorithm with 9 different ℓ 2 regularization penalty parameters (alpha): 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100. All other parameters were set to the default. • SksurvCoxPHBeamSearch (our method): 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.Recording experimental results: For each method with specific configuration, we have a set of up to 5 fitted models on each dataset. Some metrics may be unavailable:• Precision, recall, and f1-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 SksurvFastSVM 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 as corresponding error bars. Figure 7 : Optimization on the Flchain dataset with λ 1 = 1.0 and λ 2 = 1.0. The exact Newton method can be applied to solve the ℓ 1 -regularized problems, so we only compare with quasi Newton and proximal Newton. These two baselines both have the losses blow up when the ℓ 2 regularization is weak. In contrast, our methods based on the quadratic and cubic surrogate functions have losses that monotonically decrease.Figure 8 : Optimization on the Flchain dataset with λ 1 = 1.0 and λ 2 = 5.0. The exact Newton method can be applied to solve the ℓ 1 -regularized problems, so we only compare with quasi Newton and proximal Newton. Stronger ℓ 2 regularization helps these two baselines avoid the losses going into infinity. However, these two baselines are still significantly slower than our methods.Figure 9 : Optimization on the Employee Attrition dataset with λ 1 = 0 and λ 2 = 1.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 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 10 : Optimization on the Employee Attrition dataset with λ 1 = 0 and λ 2 = 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 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 : Optimization on the Employee Attrition dataset with λ 1 = 1.0 and λ 2 = 1.0. The exact Newton method can be applied to solve the ℓ 1 -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), 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.Figure 12 : Optimization on the Employee Attrition dataset with λ 1 = 1.0 and λ 2 = 5.0. The exact Newton method can be applied to solve the ℓ 1 -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), 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. The exact Newton method can be applied to solve the ℓ 1 -regularized problems, so we only compare with quasi Newton and proximal Newton. These two baselines both have the losses blow up when both the ℓ 1 and ℓ 2 regularizations are weak. In contrast, our methods based on the quadratic and cubic surrogate functions have monotonically decreasing losses.Figure 16 : Optimization on the Kickstarter1 dataset with λ 1 = 1.0 and λ 2 = 5.0. The exact Newton method can be applied to solve the ℓ 1 -regularized problems, so we only compare with quasi Newton and proximal Newton. These two baselines both have the losses blow up even when we have a stronger ℓ 2 regularization. In contrast, our methods based on the quadratic and cubic surrogate functions have monotonically decreasing losses.Figure 17 : Optimization on the Dialysis dataset with λ 1 = 0 and λ 2 = 1.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 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 18 : Optimization on the Dialysis dataset with λ 1 = 0 and λ 2 = 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 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 19 : Optimization on the Dialysis dataset with λ 1 = 1.0 and λ 2 = 1.0. The exact Newton method can be applied to solve the ℓ 1 -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), 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.Figure 20 : Optimization on the Dialysis dataset with λ 1 = 1.0 and λ 2 = 5.0. The exact Newton method can be applied to solve the ℓ 1 -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), 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. preferably report a 2-sigma error bar than state that they have a 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 of range (e.g. negative error 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.See the skglm tutorial at https://contrib.scikit-learn.org/skglm/tutorials/cox_datafit. html#maths-cox-datafit.For the ElasticNet problem where the penalty is λ1∥•∥1 + λ2∥•∥ 2 2 , 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 a λ1∥•∥ penalty. 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.
- 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
- A discussion of the tightness of the bound seems necessary when optimizing the bound instead of the loss function
questions: 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.
- 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)?
- 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: limitations: The authors have addressed some of the limitations of the model.
paper_summary: paper_summary: This paper presents a new optimization algorithm for the Cox model, which is a classical algorithm for survival analysis presented in 1972.
strength_weakness: strength_weakness: 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.

Misc:
- Line 52: we are -> ours are
- 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.
- Line 116: result result
- Line 124: a comma is missing
- Line 195: wise -> wide
- Line 257: it is better to clarify which correlation coefficient is used: Pearson, Spearman, or any other one?
- Page 15: nothing is changed in the statement "move the partial derivative operator"
- Lines 539-540: Redundant vertical lines in the inequality.
questions: questions: N/A
limitations: limitations: Discussed in page 9.
paper_summary: paper_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.
strength_weakness: strength_weakness: Excellent paper, only a few minor typos and questions, see below.
questions: questions: - It would be beneficial to add a pointer on the proof of convexity (line 82) of the nll.
- Doesn't the cumulative reverse sum require times to be sorted, i.e. an initial O(n log(n)) computation?
- Typo Line 195 - 'wise' should be 'wide'.

- 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: limitations: As an optimization paper, the authors mention that the work's limitations are the same as those of the standard Cox model.
paper_summary: paper_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.
strength_weakness: strength_weakness: I do not believe the paper makes a significant contribution to the literature.

The 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.

The 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.

I 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

I ran the following code based on the sksurv documentation:

```
from sksurv.datasets import load_flchain
from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.preprocessing import OneHotEncoder
from sklearn.impute import SimpleImputer

X, y = load_flchain() # Load the data
Xt = OneHotEncoder().fit_transform(X) # One-hot encode categorical variables
imputer = SimpleImputer().fit(Xt) # Dataset contains missing values, impute with mean
Xt = imputer.transform(Xt)

cph = CoxPHSurvivalAnalysis(alpha=1.0, verbose=1) # Fit the model with regularization lambda_2 = 1 as in Fig. 1
cph.fit(Xt, y)
```

The 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.

In addition to these major issues, there are several other problems:
1. 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.

2. 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].

3. 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).

4. 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).

5. 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".

**References:**

[1] Schoenfeld, David. Partial Residuals for the Proportional Hazards Regression Model. Biometrika (1982).
questions: questions: Is there a mathematical justification for why we should expect the method to be robust to highly correlated features?

Please also address the concerns listed in the Weaknesses section.
limitations: limitations: Limitations are briefly discussed (lines 278-282), mostly regarding limitations of the Cox model itself.

Comments: title: Follow-up on what is the benefit of our preprocessing
Maintaining score
comment: I have read the additional comments and maintain my score.
Please see the results below for the comparison between Cox models trained on the original and preprocessed features, respectively.
The results correspond to our Figure 3 and Figure 4 in the main paper.
The trend is very clear: the models trained on the preprocessed features achieve better metrics than the ones trained on the original features (for CIndex, the higher the better; for IBS, the lower the better).
Note that the Employee Attrition dataset has 17 original features while the Dialysis dataset has only 7 original features.

1. **Results on Employee Attrition; see Figure 3**

| Use Thresholding as Preprocessing | Support Size | Method | Train CIndex | Test CIndex | Train IBS | Test IBS |
| :---: | :----: | :---: | :---: | :---: | :---: | :---: |
| no | 5 | ABESS | 0.855 $\pm$ 0.002 | 0.855 $\pm$ 0.004 | 0.259 $\pm$ 0.003 | 0.260 $\pm$ 0.025 |
| yes | 5 | ABESS | 0.919 $\pm$ 0.006 | 0.918 $\pm$ 0.011 | 0.204 $\pm$ 0.010 | 0.206 $\pm$ 0.030 |
| no | 10 | ABESS | 0.856 $\pm$ 0.001 | 0.855 $\pm$ 0.007 | 0.253 $\pm$ 0.003 | 0.254 $\pm$ 0.021 |
| yes | 10 | ABESS | 0.927 $\pm$ 0.004 | 0.925 $\pm$ 0.013 | 0.199 $\pm$ 0.012 | 0.201 $\pm$ 0.032 |
| no | 15 | ABESS | 0.857 $\pm$ 0.002 | 0.856 $\pm$ 0.007 | 0.252 $\pm$ 0.002 | 0.252 $\pm$ 0.022 |
| yes | 15 | ABESS | 0.925 $\pm$ 0.007 | 0.924 $\pm$ 0.013 | 0.200 $\pm$ 0.018 | 0.201 $\pm$ 0.030 |
| no | 17 | ABESS | 0.856 $\pm$ 0.002 | 0.856 $\pm$ 0.006 | 0.252 $\pm$ 0.002 | 0.252 $\pm$ 0.022 |
| yes | 17 | ABESS | 0.924 $\pm$ 0.006 | 0.922 $\pm$ 0.011 | 0.209 $\pm$ 0.016 | 0.209 $\pm$ 0.023 |
| yes | 20 | ABESS | 0.931 $\pm$ 0.004 | 0.928 $\pm$ 0.012 | 0.189 $\pm$ 0.021 | 0.194 $\pm$ 0.039 |
| yes | 25 | ABESS | 0.935 $\pm$ 0.006 | 0.934 $\pm$ 0.007 | 0.182 $\pm$ 0.021 | 0.188 $\pm$ 0.022 |
| yes | 30 | ABESS | 0.939 $\pm$ 0.004 | 0.937 $\pm$ 0.010 | 0.172 $\pm$ 0.017 | 0.176 $\pm$ 0.034 |
| no | 5 | ours | 0.854 $\pm$ 0.002 | 0.854 $\pm$ 0.005 | 0.252 $\pm$ 0.002 | 0.252 $\pm$ 0.020 |
| yes | 5 | ours | 0.920 $\pm$ 0.001 | 0.919 $\pm$ 0.009 | 0.199 $\pm$ 0.006 | 0.201 $\pm$ 0.026 |
| no | 10 | ours | 0.854 $\pm$ 0.002 | 0.853 $\pm$ 0.007 | 0.246 $\pm$ 0.002 | 0.246 $\pm$ 0.020 |
| yes | 10 | ours | 0.943 $\pm$ 0.002 | 0.942 $\pm$ 0.007 | 0.127 $\pm$ 0.002 | 0.128 $\pm$ 0.013 |
| no | 15 | ours | 0.855 $\pm$ 0.002 | 0.854 $\pm$ 0.007 | 0.245 $\pm$ 0.002 | 0.245 $\pm$ 0.020 |
| yes | 15 | ours | 0.951 $\pm$ 0.002 | 0.950 $\pm$ 0.007 | 0.107 $\pm$ 0.003 | 0.109 $\pm$ 0.010 |
| no | 17 | ours | 0.855 $\pm$ 0.002 | 0.854 $\pm$ 0.007 | 0.245 $\pm$ 0.002 | 0.245 $\pm$ 0.020 |
| yes | 17 | ours | 0.951 $\pm$ 0.002 | 0.950 $\pm$ 0.007 | 0.105 $\pm$ 0.003 | 0.107 $\pm$ 0.010 |
| yes | 20 | ours | 0.952 $\pm$ 0.001 | 0.950 $\pm$ 0.008 | 0.103 $\pm$ 0.002 | 0.105 $\pm$ 0.009 |
| yes | 25 | ours | 0.953 $\pm$ 0.002 | 0.951 $\pm$ 0.007 | 0.101 $\pm$ 0.002 | 0.103 $\pm$ 0.007 |
| yes | 30 | ours | 0.954 $\pm$ 0.002 | 0.952 $\pm$ 0.007 | 0.099 $\pm$ 0.002 | 0.102 $\pm$ 0.007 |

2. **Results on Dialysis; see Figure 4**

| Use Thresholding as Preprocessing | Support Size | Method | Train CIndex | Test CIndex | Train IBS | Test IBS |
| :---: | :----: | :---: | :---: | :---: | :---: | :---: |
| no | 5 | ABESS | 0.649 $\pm$ 0.002 | 0.642 $\pm$ 0.008 | 0.173 $\pm$ 0.001 | 0.176 $\pm$ 0.010 |
| yes | 5 | ABESS | 0.666 $\pm$ 0.011 | 0.650 $\pm$ 0.020 | 0.178 $\pm$ 0.002 | 0.182 $\pm$ 0.012 |
| no | 7 | ABESS | 0.651 $\pm$ 0.002 | 0.646 $\pm$ 0.008 | 0.173 $\pm$ 0.001 | 0.176 $\pm$ 0.010 |
| yes | 7 | ABESS | 0.680 $\pm$ 0.008 | 0.670 $\pm$ 0.011 | 0.175 $\pm$ 0.001 | 0.178 $\pm$ 0.010 |
| yes | 10 | ABESS | 0.692 $\pm$ 0.004 | 0.683 $\pm$ 0.012 | 0.170 $\pm$ 0.002 | 0.174 $\pm$ 0.009 |
| yes | 15 | ABESS | 0.701 $\pm$ 0.008 | 0.693 $\pm$ 0.018 | 0.166 $\pm$ 0.008 | 0.170 $\pm$ 0.013 |
| yes | 20 | ABESS | 0.709 $\pm$ 0.008 | 0.694 $\pm$ 0.015 | 0.164 $\pm$ 0.003 | 0.171 $\pm$ 0.011 |
| yes | 25 | ABESS | 0.724 $\pm$ 0.013 | 0.710 $\pm$ 0.007 | 0.156 $\pm$ 0.005 | 0.163 $\pm$ 0.011 |
| yes | 30 | ABESS | 0.731 $\pm$ 0.003 | 0.719 $\pm$ 0.010 | 0.153 $\pm$ 0.003 | 0.159 $\pm$ 0.010 |
| no | 5 | ours | 0.651 $\pm$ 0.003 | 0.644 $\pm$ 0.010 | 0.179 $\pm$ 0.003 | 0.181 $\pm$ 0.010 |
| yes | 5 | ours | 0.672 $\pm$ 0.005 | 0.654 $\pm$ 0.014 | 0.176 $\pm$ 0.001 | 0.181 $\pm$ 0.001 |
| no | 7 | ours | 0.651 $\pm$ 0.003 | 0.646 $\pm$ 0.010 | 0.175 $\pm$ 0.003 | 0.177 $\pm$ 0.009 |
| yes | 7 | ours | 0.686 $\pm$ 0.003 | 0.669 $\pm$ 0.013 | 0.171 $\pm$ 0.001 | 0.176 $\pm$ 0.009 |
| yes | 10 | ours | 0.703 $\pm$ 0.007 | 0.690 $\pm$ 0.017 | 0.165 $\pm$ 0.002 | 0.170 $\pm$ 0.011 |
| yes | 15 | ours | 0.725 $\pm$ 0.011 | 0.704 $\pm$ 0.007 | 0.158 $\pm$ 0.006 | 0.164 $\pm$ 0.009 |
| yes | 20 | ours | 0.739 $\pm$ 0.005 | 0.722 $\pm$ 0.014 | 0.152 $\pm$ 0.003 | 0.157 $\pm$ 0.010 |
| yes | 25 | ours | 0.742 $\pm$ 0.004 | 0.722 $\pm$ 0.014 | 0.150 $\pm$ 0.002 | 0.156 $\pm$ 0.011 |
| yes | 30 | ours | 0.745 $\pm$ 0.003 | 0.719 $\pm$ 0.014 | 0.149 $\pm$ 0.002 | 0.157 $\pm$ 0.010 |
1. **Your claim that using Eq.8 and standard methods for reducing runtime when computing multiple partial sums, the runtime reduction is not surprising**
    
Your opinion on runtime reduction not being surprising is based on the assumption that you already have access to Eq.8.
However, before our submission, no one knows we could possibly write the 2nd-order partial derivatives in such a way.
Eq.8 is unique to the Cox model (Eq.8) and does not hold for a generic loss function.
**Since Eq.8 is not something already known/established, it's not possible to conclude our work ''is not surprising'' considering only what follows it**.
As mentioned before, the previous approaches rely on $\boldsymbol{\eta}^T \nabla_{\boldsymbol{\eta}} \ell(\boldsymbol{\eta}) \boldsymbol{\eta}$.
This is why we can reduce the complexity from $O(n^2)$ to $O(n)$ and why we say the runtime reduction is surprising.

2. **What is the benefit of our preprocessing?**

We do not perform this for the sake of creating a challenging dataset.
It does have a potential benefit for producing better models.
Due to the character limit, please see another post where we support this claim through experiments.
    
The origin of our preprocessing can be traced back to additive models in statistics (see Line 27, 257-258, and 280, with citations [1, 11, 27, 28, 49, 67]).
Additive models have been extensively used to capture nonlinear relationships between the target and continuous variables.
For example, on the ICU patients, death risks may not increase linearly w.r.t age, a continuous variable; the risk for older patients may increase at a much higher rate.
Using the original continuous variable would fail to capture this relationship.
    
Lastly, additive models can be created by other preprocessing procedures, such as splines and polynomials (see [27, 28, 67]).
However, no matter which procedure we use, the preprocessed features are likely much more correlated than the original ones.
As shown in Section 4.2 (variable selection), our method can handle highly correlated features much better than existing methods.
    
3. **Minimizing a cubic by finding the roots of a quadratic equation still does not constitute significant mathematical novelty**

We try to minimize the following cubic surrogate function:
$$h_x(\Delta x):= f(x) + f'(x) \Delta x + \frac{1}{2} f''(x) \Delta x^2 + \frac{1}{6} L_3 \vert \Delta x \vert^3.$$
Your opinion on our method not constituting significant mathematical novelty is based on the assumption that *$h_x(\Delta x)$ already exists*.
Specifically, you assume that $f'(x)$, $f''(x)$, and $L_3$ have already been given to you.
However, this is not the case, due to two things mentioned in our initial rebuttal:

a) $f''(x)$can be computed in $O(n)$ time. Previously, the best way takes $O(n^2)$. *If you compute $f''(x)$ in the old-fashioned way, it takes much longer to formulate $h_x(\Delta x)$*.

b) Previously, it is unknown, regarding the Cox model, whether $L_3$ (Lipschitz constant for $f''(x)$) exists and, if so, whether there exists an explicit way to compute $L_3$. In our work, we answer both questions positively, through Eq.9 and Eq.14. *Without our work, you cannot even compute $L_3$ and thus formulate $h_x(\Delta x)$, let alone minimize it*.

Moreover, *obtaining Eq.9 and Eq.14 are nontrivial mathematical tasks*.
No one has expected these nice structures (probabilistic formula for the 3rd-order partial derivatives and its connection to the 3rd-order central moment calculation) could pop up out of the Cox model.

Lastly, there are subjective opinions and objective facts.
Regarding subjective opinions, we accept that different researchers can have different research tastes on novelty.
However, objectively, the fact is that no one has written Eq.18, the analytical minimizer to $h_x(\Delta x)$. If you think they have, please point us to the literature where they publish this equation.

4. **Even with the addition of L1 regularization, this is reasonably behaved one-dimensional optimization problem which should not pose much challenge for standard approaches**

Our reasoning above also applies here.
*Without our work, you cannot even formulate this L1-regularized cubic surrogate function, let alone minimize it.*
Previous researchers didn't know that the cubic surrogate function could be useful to solving the L1-regularized problem.
All they had are the quasi and proximal Newton methods.
Again, we make a contribution here because Eq.22 has never been presented.
Moreover, our Eq.22 is one step further beyond the traditional soft-thresholding operator (well known as solving the L1-regularized problem with the quadratic surrogate function).
Our work opens the door to solving the L1-regularized problem with the cubic surrogate function.
Thank you for your answers, I am maintaining my score
Thanks to the authors for their detailed response. My questions related to points 3 & 4, as well as the monotonic decrease in the objective have all been answered, and I have raised my score accordingly. I still lean towards rejection for the following reasons:

1. I understand that naively computing the second order derivatives naively using the formula on line 146 would lead to a high runtime. However, using equation (8) and standard methods for reducing runtime when computing multiple partial sums, the runtime reduction is not surprising.

2. I now understand the reason for the difference in wall clock runtime from my implementation vs. yours, and I believe that the baselines were implemented correctly. However, I am not convinced that the preprocessing used in the paper is realistic. Turning continuous features into highly correlated categorical features is the opposite of what one would usually do in practice. Is there any benefit to this procedure? For instance, does it improve the downstream performance of the model, e.g. in terms of the C-index on the test data? Barring some good justification for this (other than that it leads to better performance for the proposed method over existing baselines), this experiment doesn't seem convincing.

6. I agree with the authors that the cubic equation does not apply to this case and I apologize for the misunderstanding. That said, minimizing a cubic by finding the roots of a quadratic equation still does not constitute significant mathematical novelty. Even with the addition of L1 regularization, this is reasonably behaved one-dimensional optimization problem which should not pose much challenge for standard approaches.
Thank you for your comments. I will maintain my original score.
Thank you for the response, I have increased the score to accept.