Paper deep dive
Distributionally Robust Survival Models under Subpopulation Shift and Outlier Contamination
Seonghwi Kim, Sung Ho Jo, Minwoo Chae
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Learning robust survival models under distribution shift is an important but challenging problem in many applications. In heterogeneous populations, a model that performs well on average may still perform poorly on certain subpopulations, and this issue becomes even more severe when the training data are contaminated by outliers. In this paper, we propose a novel distributionally robust framework for survival analysis that jointly addresses latent subpopulation shift and outlier contamination. The proposed method combines an outer minimization that selects a refined nominal distribution by reducing the influence of contaminated samples and an inner maximization that focuses on the most challenging subpopulation. This formulation directly accommodates non-decomposable survival losses while preserving interactions across samples, including the risk-set structure of the Cox negative partial log-likelihood. We develop an alternating gradient-based algorithm with outer updates derived from the KKT conditions of the inner maximization. Experiments on simulated data and two survival benchmarks demonstrate that the proposed method remains robust when subpopulation shift and outlier contamination occur simultaneously. It stabilizes training in contaminated settings and substantially improves worst-group performance across both linear and nonlinear survival models, while maintaining competitive and sometimes superior overall performance.
Tags
Links
- Source: https://arxiv.org/abs/2610.02868v1
- Canonical: https://arxiv.org/abs/2610.02868v1
Trouble viewing inline? Open PDF directly →
Full Text
128,453 characters extracted from source content.
Expand or collapse full text
Distributionally Robust Survival Models under Subpopulation Shift and Outlier Contamination Seonghwi Kim Affiliation: Pohang University of Science and Technology Email: kshwi@postech.ac.kr Sung Ho Jo Affiliation: Pohang University of Science and Technology Email: tjdgh1813@postech.ac.kr Minwoo Chae †thanks: Corresponding author. Affiliation: Pohang University of Science and Technology Email: mchae@postech.ac.kr Abstract Learning robust survival models under distribution shift is an important but challenging problem in many applications. In heterogeneous populations, a model that performs well on average may still perform poorly on certain subpopulations, and this issue becomes even more severe when the training data are contaminated by outliers. In this paper, we propose a novel distributionally robust framework for survival analysis that jointly addresses latent subpopulation shift and outlier contamination. The proposed method combines an outer minimization that selects a refined nominal distribution by reducing the influence of contaminated samples and an inner maximization that focuses on the most challenging subpopulation. This formulation directly accommodates non-decomposable survival losses while preserving interactions across samples, including the risk-set structure of the Cox negative partial log-likelihood. We develop an alternating gradient-based algorithm with outer updates derived from the KKT conditions of the inner maximization. Experiments on simulated data and two survival benchmarks demonstrate that the proposed method remains robust when subpopulation shift and outlier contamination occur simultaneously. It stabilizes training in contaminated settings and substantially improves worst-group performance across both linear and nonlinear survival models, while maintaining competitive and sometimes superior overall performance. 1 Introduction Survival analysis is a fundamental problem in statistics and machine learning, with broad applications in healthcare, reliability analysis, and customer retention. In many applications, the population comprises subpopulations with distinct survival patterns. Such heterogeneity may arise across demographic groups such as age, sex, and race, as well as across hospitals, clinical environments, or latent patient subtypes (Hu and Chen, 2024; Rahman and Purushotham, 2022; Fan et al., 2024). When the mixture proportions of these subpopulations differ between training and test distributions, a model with strong average performance can still perform poorly on some subpopulations (Sagawa et al., 2020; Yang et al., 2023; Hu and Chen, 2024). To address this challenge, distributionally robust learning, or distributionally robust optimization11 1 While DRO was originally studied in the context of optimization, we use the term more broadly to refer to learning methods for robustness under distribution shift. (DRO), has been widely studied as a principled framework for improving robustness under distributional shift (Mohajerin Esfahani and Kuhn, 2018; Duchi and Namkoong, 2021). Applying DRO to survival analysis, however, is not straightforward. Many commonly used survival losses are non-decomposable, involving interactions across samples. For example, the Cox negative partial log-likelihood depends on risk-set interactions, whereas ranking-based survival losses are defined through pairwise comparisons across samples (Cox, 1972; Raykar et al., 2007; Lee et al., 2018). In contrast, many existing DRO methods assume additive losses, which decompose into a sum of terms that each depend on a single sample. This structural difference prevents these methods from being directly applied to such survival losses (Jin et al., 2025; Hu and Chen, 2022; Hu and Chen, 2024). Recently, Hu and Chen (2024) proposed an f-divergence based DRO (f-DRO) framework for survival models, focusing on the χ2χ^2-divergence, and introduced a sample-splitting strategy to handle the interaction structure of survival losses. Their results show that f-DRO can improve robustness to latent subpopulation shift in survival analysis. However, this approach is highly sensitive to outliers or contaminated observations, since it improves robustness by assigning larger weights to high-loss samples. DORO (Zhai et al., 2021) addresses both subpopulation shift and outlier contamination, but is designed for additive losses and does not directly apply to non-decomposable survival losses with cross-sample dependence. To address these limitations, we propose a coupled distributionally robust framework for survival analysis that jointly addresses latent subpopulation shift and outlier contamination. The proposed formulation is based on a mixture-based ambiguity set with different roles in a coupled minimax problem. The outer minimization refines the nominal distribution to reduce the influence of potentially contaminated samples. Based on this refined distribution, the inner maximization concentrates on the most challenging subpopulation. In this way, the proposed method iteratively refines the nominal distribution by reducing the influence of contaminated samples and reweights the most challenging subpopulation within the refined distribution, providing robustness to both outlier contamination and subpopulation shift. For differentiable non-decomposable survival losses, we develop an efficient alternating gradient-based algorithm with outer updates derived from the KKT conditions of the inner maximization. We instantiate the framework with the Cox negative partial log-likelihood, preserving its risk-set interactions. We further show that, for additive losses, the formulation admits a two-stage structure: the outer minimization performs ϵε-trimming, and the inner maximization reweights the highest-loss α-fraction of the remaining probability mass. Through simulation studies and experiments on two widely used survival benchmarks, METABRIC and FLC, we show that the proposed method remains robust under subpopulation shift and outlier contamination. In simulations, it exhibits stable training trajectories even in contaminated settings, whereas f-DRO exhibits unstable trajectories. On real benchmarks, it yields substantially improved worst-group performance across both linear and nonlinear survival models. It also achieves competitive overall performance, often even outperforming existing baselines in contaminated settings. Our main contributions can be summarized as follows. • We propose a novel coupled distributionally robust framework for survival analysis that jointly addresses two practically important challenges: latent subpopulation shift and outlier contamination. • We develop an alternating gradient-based algorithm that directly optimizes non-decomposable survival losses while preserving interactions across samples. Its outer updates are derived from the KKT conditions of the inner maximization. • We demonstrate through simulation studies and real data experiments that the proposed method remains stable under outlier contamination and achieves substantially improved worst-group performance under subpopulation shift. 2 Related Work Survival Analysis under Subpopulation Shift Survival analysis has traditionally been developed based on classical statistical frameworks, including nonparametric estimators, proportional hazards models, and related semiparametric approaches (Cox, 1972; Kaplan and Meier, 1958; Buckley and James, 1979; Bennett, 1983; Kalbfleisch and Prentice, 2002). More recently, it has been extended to machine learning methods that allow more flexible nonlinear modeling of survival outcomes (Ishwaran et al., 2008), and further to deep learning based approaches designed for richer representation learning and more flexible risk prediction (Lee et al., 2018; Ranganath et al., 2016; Katzman et al., 2018; Avati et al., 2020; Tang et al., 2023). However, these models do not explicitly account for subpopulation shift, and therefore often perform poorly on specific groups. This issue has been studied in survival analysis mainly from two directions: fairness across observed groups and federated survival analysis under institutional heterogeneity. Fairness-based approaches typically impose constraints or regularization to encourage similar predictions or predictive performance across demographic or clinically defined groups (Keya et al., 2021; Zhang and Weiss, 2022; Zhang and Weiss, 2023; Rahman and Purushotham, 2022; Do et al., 2023). Federated survival methods, in contrast, address heterogeneity across hospitals or institutions with different local patient populations and data distributions (Andreux et al., 2020; Archetti et al., 2023; Seidi et al., 2024; Archetti and Matteucci, 2023; Rahman and Purushotham, 2023; Rahman and Purushotham, 2025; Wang et al., 2022). However, these approaches generally rely on explicit group information, institution-level partitioning, or other prespecified heterogeneous structure. Distributionally Robust Optimization Distributionally robust learning has emerged as a promising framework for handling various forms of distributional shift (Lin et al., 2022; Delage and Ye, 2010; Mohajerin Esfahani and Kuhn, 2018; Kuhn et al., 2019; Gao et al., 2024; Sagawa et al., 2020). In particular, DRO has been widely studied for latent subpopulation shift, where robustness must be achieved without explicit group annotations. A prominent line of work improves worst-case performance by emphasizing high-loss regions of the distribution, as in f-DRO and related CVaR-type objectives (Namkoong and Duchi, 2016; Duchi and Namkoong, 2021; Hashimoto et al., 2018; Shapiro et al., 2021). Recent work has extended this idea to survival analysis without relying directly on predefined sensitive attributes (Hu and Chen, 2022; Hu and Chen, 2024). However, these methods can be highly vulnerable to outliers and may become unstable under contamination (Zhai et al., 2021). This issue is particularly important in survival analysis, where outliers are pervasive and can substantially affect model fitting and prediction. Existing robust survival methods mitigate outlier effects (Heritier et al., 2009; Sinha, 2019; Farcomeni and Ventura, 2012; Luo et al., 2022; Farcomeni and Viviani, 2011), but do not directly address latent subpopulation shift. Zhai et al. (2021) improve outlier robustness in f-DRO but focus on additive losses and do not cover non-decomposable survival losses with cross-sample dependence. Thus, jointly addressing outlier contamination and latent subpopulation shift remains an important problem in survival analysis. Independently of our work, Huang et al. (2026) propose a nested DRO framework for robust learning with additive losses. Their formulation does not directly apply to non-decomposable survival losses with cross-sample dependence. Other Approaches to Subpopulation Shift Subpopulation shift has also been studied through latent subgroup discovery and annotation-free group robustness (Liu et al., 2021; Hong et al., 2024; Stromberg et al., 2024), group-aware priors (Rudner et al., 2024), outlier-robust group inference (Zeng et al., 2023), and DRO frameworks with hierarchical or structured ambiguity sets (Jo et al., 2026; Kim et al., 2026). Our work focuses on jointly addressing subpopulation shift and outlier contamination in survival models with non-decomposable losses. 3 Problem Setup In this section, we introduce the notation and formalize the problem setup. Although our primary interest lies in survival analysis, we first present our framework in a general statistical learning setting, and then discuss its application to survival models in Section 4.3. Let Z∈Z denote a random variable. For any distribution Q, we denote expectation under Q by Q[⋅]E_Q[·]. For a positive integer n, let [n]:=1,…,n[n]:=\1,…,n\ denote the index set, and let n−1:=p∈ℝn:pi≥0,∑i=1npi=1S_n-1:=\p ^n:p_i≥ 0,\ _i=1^np_i=1\ denote the (n−1)(n-1)-simplex. Given a model fθf_θ parameterized by θ, we denote the corresponding loss function by ℓ(Z,θ) (Z,θ). We consider a learning problem under subpopulation shift with outlier contamination. Specifically, the target distribution PtestP_ test is assumed to be a mixture of K latent groups, Ptest=∑k=1KπkPkP_ test= _k=1^K _kP_k, where ∑k=1Kπk=1 _k=1^K _k=1 and πk∈(0,1) _k∈(0,1). Here, PkP_k denotes the distribution of the k-th latent group, and πk _k denotes its mixture proportion. Let PtrainP_train denote the clean training distribution before outlier contamination, which may have different mixing proportions from PtestP_ test. The observed training distribution follows Huber’s ϵε-contamination model (Huber, 1992), P~train=(1−ϵ)Ptrain+ϵP~ P_train=(1-ε)P_train+ε P, where P~ P is an arbitrary outlier distribution and ϵ∈[0,1)ε∈[0,1) denotes the contamination level. Let =zii=1nD=\z_i\_i=1^n be a dataset drawn i.i.d. from P~train P_train. Our goal is to learn a model that is robust to both outlier contamination in the training data and subpopulation shift under the latent mixture distribution PtestP_ test. A classical baseline is empirical risk minimization (ERM), which minimizes the average loss over the observed training data. However, ERM does not account for distributional shift between the clean training distribution PtrainP_train and the target distribution PtestP_ test, or for outlier contamination in the observed training distribution P~train P_train. A promising alternative is distributionally robust learning, which minimizes the worst-case risk over a set of plausible distributions. Specifically, it solves minθmaxQ∈ℬρ(P¯)Q[ℓ(Z,θ)], _θ _Q _ρ( P)\;E_Q [ (Z,θ) ], where P¯ P is a nominal distribution and ℬρ(P¯)B_ρ( P) is an ambiguity set that captures distributional uncertainty around P¯ P. Two prominent DRO formulations for subpopulation shift are Group DRO (Sagawa et al., 2020) and f-DRO (Namkoong and Duchi, 2016; Duchi and Namkoong, 2021). Group DRO directly optimizes the worst-case group risk but requires access to group labels, which are unavailable in many practical settings, including ours. In contrast, f-DRO defines the ambiguity set via an f-divergence and achieves robustness by reweighting high-loss samples, which can effectively address latent subpopulation shift. However, it remains inherently vulnerable to outliers (Zhai et al., 2021). Therefore, it may fail to provide robustness when both subpopulation shift and outlier contamination are present. 4 Proposed Method We first introduce the general coupled DRO formulation and illustrate its structure with additive losses. We then extend the framework to non-decomposable survival losses, the main focus of this work. 4.1 Formulation To motivate the proposed formulation, we first define a mixture ambiguity set. For a given nominal distribution P¯ P and a constant γ∈(0,1]γ∈(0,1], define γ(P¯)=Q:P¯=γQ+(1−γ)Q~for someQ~∈(),Q_γ( P)= \Q: P=γ Q+(1-γ) Q\; for some\; Q (Z) \, (1) where ()P(Z) denotes the set of all probability distributions on Z. This set consists of γ-mixture components of P¯ P, each of which is absolutely continuous with respect to P¯ P, denoted by Q≪P¯Q P. This ambiguity set serves different roles in minimization and maximization. For a fixed θ, consider minQ∈1−ϵ(P¯)Q[ℓ(Z,θ)]andmaxQ∈α(P¯)Q[ℓ(Z,θ)]. _Q _1-ε( P)E_Q[ (Z,θ)] _Q _α( P)E_Q[ (Z,θ)]. (2) Here, ϵ∈[0,1)ε∈[0,1) and α∈(0,1]α∈(0,1]. The minimization selects a distribution Q∈1−ϵ(P¯)Q _1-ε( P) that minimizes the expected loss by assigning smaller weights to high-loss regions, including potential outliers. The resulting distribution can be viewed as a refined version of P¯ P, with ϵε reflecting the degree of potential outlier contamination. In contrast, the maximization selects a distribution Q∈α(P¯)Q _α( P) that maximizes the expected loss by assigning larger weights to high-loss regions. This distribution emphasizes regions that are particularly difficult under the current model, with α controlling the size of the subpopulation considered. This maximization is closely related to CVaR-type DRO (Shapiro et al., 2021), which focuses on the high-loss tail of the distribution. Motivated by these two roles, we propose the following coupled DRO problem: minθ[minQ′∈1−ϵ(P¯)maxQ∈α(Q′)Q[ℓ(Z,θ)]]. _θ [ _Q _1-ε( P) _Q _α(Q )E_Q [ (Z,θ) ] ]. (3) We refer to the minimization over Q′Q and the maximization over Q as the outer minimization and the inner maximization, respectively. The formulation couples nominal refinement with adversarial reweighting through the dependence of the inner ambiguity set α(Q′)Q_α(Q ) on the refined distribution Q′Q . 4.2 Structural analysis under additive losses To facilitate understanding of the proposed formulation, we consider a class of losses for which the optimization structure becomes explicit. Specifically, given a weight vector q∈n−1q _n-1 over the dataset D, we view the empirical risk as a function ℒ:n−1×Θ→ℝL:S_n-1× . We call this class of losses additive if the empirical risk can be decomposed as ℒ(q,θ)=∑i=1nqiℓ(zi,θ)L(q,θ)= _i=1^nq_i (z_i,θ), where each loss term depends on a single sample. This class includes standard supervised learning losses such as the cross-entropy loss and the squared loss. We now take the nominal distribution P¯ P to be supported on D with probability mass vector p¯∈n−1 p _n-1. Under this choice, each feasible distribution can be represented by its probability mass vector over D. Thus, for additive losses, the proposed formulation in (3) admits the following equivalent discrete representation: minθ[minq′∈1−ϵ(p¯)maxq∈α(q′)ℒ(q,θ)]. _θ [ _q _1-ε( p) _q _α(q )L(q,θ) ]. (4) Here, for a reference weight vector p∈n−1p _n-1 and γ∈(0,1]γ∈(0,1], the ambiguity set γ(p)Q_γ(p) is defined as γ(p)=q∈n−1:0≤qi≤piγ,∀i∈[n].Q_γ(p)= \q _n-1:0≤ q_i≤ p_iγ,\;∀ i∈[n] \. (5) Details of this finite-support representation of the ambiguity set are provided in Lemma 1. The discrete coupled formulation in (4) admits an explicit two-stage characterization. Specifically, the outer minimization performs ϵε-trimming, and the inner maximization reweights the highest-loss α-fraction of the remaining probability mass. A precise characterization and the corresponding procedure are given in Appendix A.1.1 and Algorithm 2. Thus, our formulation includes CVaR-DORO (Zhai et al., 2021) as a special case for additive losses. However, this structural result relies critically on the additivity of the loss and does not extend to settings where the loss depends on multiple samples simultaneously. In this paper, we extend the proposed framework beyond the additive setting to handle such non-decomposable losses. Survival analysis provides a primary example of this class, which we study in the following section. 4.3 Non-decomposable survival losses In this section, we focus on right-censored survival data and specify the corresponding survival loss, although the framework can be adapted to other censoring or truncation mechanisms. Let Z=(X,Y,Δ)Z=(X,Y, ) be a random variable, where X∈X denotes the covariate, T is the survival time, and C is the censoring time. The observed time is given by Y=minT,CY= \T,C\, and the event indicator is Δ=(T≤C) =I(T≤ C). Each observed sample is denoted by zi=(xi,yi,δi)z_i=(x_i,y_i, _i). Survival losses such as the Cox negative partial log-likelihood and ranking-based losses involve cross-sample dependence. To apply the proposed coupled DRO formulation, we express the empirical survival risk as ℒ(q,θ)L(q,θ) to make its dependence on the sample weights explicit. As a representative example, we consider the Cox proportional hazards model (Cox, 1972). Given q∈n−1q _n-1, the empirical risk is defined as the weighted negative partial log-likelihood: ℒ(q,θ)=−∑i=1nqiδi[fθ(xi)−log(∑j∈[n]:yj≥yiqjexp(fθ(xj)))].L(q,θ)=- _i=1^nq_i _i [f_θ(x_i)- ( _j∈[n]:y_j≥ y_iq_j (f_θ(x_j)) ) ]. We incorporate this loss into the coupled formulation in (4). Because the sample weights also appear inside the logarithmic risk-set terms, ℒ(q,θ)L(q,θ) is nonlinear in q. Thus, the closed-form trimming and reweighting structure derived for additive losses does not carry over directly. Moreover, standard DRO formulations based on dual reformulations are difficult to apply directly to non-decomposable survival losses, since they typically rely on additive losses (Hu and Chen, 2024; Jin et al., 2025). We therefore develop a gradient-based optimization scheme that preserves the original survival loss structure with cross-sample dependence. 5 Algorithm In this section, we present the learning algorithm and its intuition. Throughout this section, we assume that ℒ(q,θ)L(q,θ) is continuously differentiable in both q and θ, which enables gradient-based updates. Our approach treats the optimization as a three-block problem over (q,q′,θ)(q,q ,θ) and adopts an alternating update strategy. A key challenge lies in updating q′q , because the outer objective is defined implicitly through the inner maximization Φα(q′,θ):=maxq∈α(q′)ℒ(q,θ) _α(q ;θ):= _q _α(q )L(q,θ). Thus, its gradient with respect to q′q is not directly available. The following proposition characterizes a supergradient through the KKT multipliers of the inner problem. Proposition 1 (Outer supergradient via KKT multipliers). Fix θ and let q⋆q be a global maximizer of the inner problem. Under the regularity and local stability assumptions stated in Appendix A.1.2, any KKT multiplier vector μ⋆μ for qi≤qi′/αq_i≤ q _i/α at q⋆q yields a supergradient μ⋆/αμ /α of Φα(⋅,θ) _α(·;θ) at q′q : Φα(q′+h,θ)≤Φα(q′,θ)+⟨μ⋆α,h⟩+o(‖h‖). _α(q +h;θ)≤ _α(q ;θ)+ μ α,h +o(\|h\|). Here, q′+h∈n−1q +h _n-1 and o(‖h‖)/‖h‖→0o(\|h\|)/\|h\|→ 0 as h→0h→ 0. Based on Proposition 1, we can derive a supergradient with respect to q′q from the KKT multipliers μ⋆μ of the inner problem. The detailed characterization of these multipliers and the derivation of this result are provided in Appendix A.1.2. Combining the q′q update based on this supergradient with gradient steps for q and θ yields the alternating procedure in Algorithm 1. We detail the full-batch updates below for clarity; stochastic mini-batch implementations are also possible (Kvamme et al., 2019; Zeng et al., 2026). Algorithm 1 Learning algorithm for the coupled DRO objective 1: Input: dataset D; step sizes ηq,ηq′,ηθ _q, _q , _θ; update steps T, K; initial values θ(0),q(0),q′(0)θ^(0),q^(0),q (0) 2: for k=0,1,2,…,K−1k=0,1,2,…,K-1 do 3: Set q(k,0)←q(k)q^(k,0)← q^(k), q′(k,0)←q′(k)q (k,0)← q (k) ⊳ Update (q,q′)(q,q ) 4: for t=0,…,T−1t=0,…,T-1 do 5: Update q(k,t+1)q^(k,t+1) using the exponentiated gradient ascent and projection in (6) 6: Update q′(k,t+1)q (k,t+1) using the exponentiated gradient descent and projection in (7) 7: end for 8: Set q(k+1)←q(k,T)q^(k+1)← q^(k,T), q′(k+1)←q′(k,T)q (k+1)← q (k,T) ⊳ Update θ 9: Update θ(k+1)θ^(k+1) using the gradient descent step in (8) 10: end for Update of q. Given θ and q′q , we update q to q q using the projected exponentiated gradient ascent update (Sagawa et al., 2020) q~i q_i ←qiexp(ηq[∇qℒ(q,θ)]i)∑j=1nqjexp(ηq[∇qℒ(q,θ)]j),i=1,…,n, ← q_i ( _q[ _qL(q,θ)]_i ) _j=1^nq_j ( _q[ _qL(q,θ)]_j ), i=1,…,n, (6) q q ←Πα(q′)(q~), ← _Q_α(q )( q), where ∇qℒ(q,θ) _qL(q,θ) denotes the gradient of the map q↦ℒ(q,θ)q (q,θ), ηq>0 _q>0 is a step size, Π(⋅) _A(·) denotes projection onto the set A, and [x]j[x]_j denotes the j-th component of a vector x. Update of q′q . Given θ and q q, we update q′q to q^′ q using projected exponentiated descent with a supergradient estimate s: q~i′ q _i ←qi′exp(−ηq′si)∑j=1nqj′exp(−ηq′sj),i=1,…,n, ← q _i (- _q s_i ) _j=1^nq _j (- _q s_j ), i=1,…,n, (7) q^′ q ←Π1−ϵ(p¯)(q~′), ← _Q_1-ε( p)( q ), where ηq′>0 _q >0 is a step size. For a nonempty interior set ℐI, we compute s from q q as follows: ℐ:=i∈[n]:0<q^i<1αqi′,ν:=1|ℐ|∑i∈ℐ[∇qℒ(q^,θ)]i,I:= \i∈[n]:0< q_i< 1αq _i \, ν:= 1|I| _i [ _qL( q,θ)]_i, si=1α([∇qℒ(q^,θ)]i−ν)+,if 0<q^i=1αqi′,0,otherwise.s_i= cases 1α ([ _qL( q,θ)]_i-ν )_+,&if 0< q_i= 1αq _i,\\[4.0pt] 0,&otherwise. cases If ℐ=∅I= , we leave q′q unchanged, as in Algorithm 3. Figure 1: Test C-index trajectories under subpopulation shift with outlier contamination. (Left) Overall test C-index for the Cox model, f-DRO, and the proposed method with (α,ϵ)=(0.1,0)(α,ε)=(0.1,0). (Middle) Overall test C-index for the proposed method with (α,ϵ)=(0.1,0.05)(α,ε)=(0.1,0.05). (Right) Worst-group test C-index for the proposed method with (α,ϵ)=(0.1,0.05)(α,ε)=(0.1,0.05). For the proposed method, only ϵε is varied from 00 to 0.050.05, while all other training configurations are kept fixed. Update of θ. After alternating the updates of q and q′q for several iterations, we update θ to θ θ using the resulting q q via a gradient descent step. θ^←θ−ηθ∇θℒ(q^,θ), θ←θ- _θ _θL( q,θ), (8) where ηθ>0 _θ>0 is a step size. To interpret the overall procedure, recall the finite-support representation of our ambiguity set: 1−ϵ(p¯) _1-ε( p) =q′∈n−1:0≤qi′≤p¯i1−ϵ,∀i∈[n], = \q _n-1:0≤ q _i≤ p_i1-ε,\;∀ i∈[n] \, α(q′) _α(q ) =q∈n−1:0≤qi≤qi′α,∀i∈[n]. = \q _n-1:0≤ q_i≤ q _iα,\;∀ i∈[n] \. The bound qi′/αq_i /α determines the maximum weight that can be assigned to sample i. Within these bounds, the inner update uses the gradient of the full survival loss to increase the loss. The outer update uses KKT multipliers from the inner problem to update q′q , thereby adjusting the upper bounds qi′/αq_i /α on the sample weights. Both updates account for how changes in sample weights affect the full loss, including other samples’ risk-set terms in the Cox model. Thus, q and q′q can be updated while preserving the cross-sample interactions in the survival loss. Computational complexity is analyzed in Appendix A.1.3. A local convergence analysis is provided in Appendix A.1.4. 6 Experiments 6.1 Simulation We evaluate robustness under subpopulation shift with outlier contamination; results without outliers are provided in Appendix A.3.1. For f-DRO, we use the χ2χ^2-divergence and sample splitting following Hu and Chen (2024). The training distribution mixes two groups with an outlier component: P~train=π1trP1+π2trP2+π~trP~ P_ train=π^tr_1P_1+π^tr_2P_2+ π^tr P, with (π1tr,π2tr,π~tr)=(0.9,0.08,0.02)(π^tr_1,π^tr_2, π^tr)=(0.9,0.08,0.02) and Ntrain=3,000N_train=3,000. The test distribution is clean and is given by Ptest=π1teP1+π2teP2P_ test=π^te_1P_1+π^te_2P_2, with (π1te,π2te)=(0.8,0.2)(π^te_1,π^te_2)=(0.8,0.2) and Ntest=5,000N_test=5,000. For each group g∈1,2g∈\1,2\, covariates are generated as X∼(0,2I2)X (0,2I_2), and survival times follow an exponential hazard model hg(t∣x)=h0exp(x⊤βg)h_g(t x)=h_0 (x _g). Censoring times are drawn from an exponential distribution, yielding 10% censoring. We set β1=(1,1) _1=(1,1) and β2=(1,−1) _2=(1,-1) for the two groups, and use β~=(−2,−2) β=(-2,-2) for the outlier component. Figure 1 compares overall and worst-group test C-index trajectories, where the latter is the minimum across the two groups (Sagawa et al., 2020). Both f-DRO and the proposed method with (α,ϵ)=(0.1,0)(α,ε)=(0.1,0) exhibit large fluctuations, consistent with their emphasis on high-loss observations under contamination. The standard Cox model remains relatively stable but has a worst-group C-index of approximately 0.550.55. Setting ϵ=0.05ε=0.05 stabilizes the proposed method’s overall and worst-group trajectories and yields a worst-group C-index near 0.700.70, illustrating the benefit of controlling outlier influence while protecting minority-group performance. Further simulation analyses, including nonlinear risk functions, higher-dimensional covariates, smaller training samples, heavier censoring, and IBS trajectories, are provided in Appendices A.3.2–A.3.4. Sensitivity to the robustness parameters α and ϵε is examined in Appendix A.2.1. Table 1: Results on METABRIC under subpopulation shift with outlier contamination. The 10% and 5% settings specify the validation tail loss for Robust-Cox and the DRO-based methods; Elastic-Cox and Cox-N use validation log-likelihood. Best results are in bold, and second-best results are underlined. Values are mean (standard deviation) over 10 independent runs. Model Method IBS ↓ Test log-likelihood ↑ Overall C-index ↑ Worst-group C-index ↑ Test tail loss ↓ Linear (10%) Elastic-Cox 0.176 (0.006) -5.061 (0.035) 0.618 (0.022) 0.501 (0.058) 6.305 (0.111) Robust-Cox 0.167 (0.005) -5.016 (0.028) 0.665 (0.017) 0.562 (0.042) 6.229 (0.048) f-DRO 0.183 (0.005) -5.121 (0.033) 0.667 (0.015) 0.574 (0.030) 5.825 (0.015) Exact-Cox 0.184 (0.006) -5.121 (0.035) 0.583 (0.028) 0.459 (0.057) 6.179 (0.164) Ours 0.184 (0.005) -5.125 (0.034) 0.684 (0.015) 0.613 (0.028) 5.835 (0.011) Linear (5%) Elastic-Cox 0.176 (0.006) -5.061 (0.035) 0.618 (0.022) 0.501 (0.058) 6.305 (0.111) Robust-Cox 0.168 (0.006) -5.020 (0.026) 0.658 (0.022) 0.547 (0.049) 6.240 (0.056) f-DRO 0.186 (0.003) -5.141 (0.018) 0.652 (0.033) 0.558 (0.046) 5.824 (0.014) Exact-Cox 0.184 (0.006) -5.121 (0.035) 0.583 (0.028) 0.459 (0.057) 6.179 (0.164) Ours 0.187 (0.004) -5.145 (0.026) 0.680 (0.015) 0.605 (0.030) 5.831 (0.007) Nonlinear (10%) Cox-N 0.181 (0.006) -5.101 (0.050) 0.601 (0.030) 0.508 (0.091) 6.561 (0.310) Robust-Cox 0.175 (0.005) -5.063 (0.032) 0.635 (0.028) 0.545 (0.049) 6.295 (0.161) f-DRO 0.182 (0.004) -5.119 (0.027) 0.644 (0.036) 0.533 (0.047) 5.839 (0.041) Exact-Cox 0.197 (0.007) -5.210 (0.049) 0.456 (0.091) 0.395 (0.078) 6.226 (0.178) Ours 0.176 (0.004) -5.074 (0.027) 0.661 (0.020) 0.585 (0.040) 5.902 (0.014) Nonlinear (5%) Cox-N 0.181 (0.006) -5.101 (0.050) 0.601 (0.030) 0.508 (0.091) 6.561 (0.310) Robust-Cox 0.175 (0.005) -5.063 (0.032) 0.635 (0.028) 0.545 (0.049) 6.295 (0.161) f-DRO 0.186 (0.004) -5.147 (0.029) 0.627 (0.036) 0.511 (0.051) 5.834 (0.021) Exact-Cox 0.197 (0.007) -5.210 (0.049) 0.456 (0.091) 0.395 (0.078) 6.226 (0.178) Ours 0.177 (0.005) -5.079 (0.033) 0.656 (0.017) 0.570 (0.036) 5.901 (0.026) 6.2 Real data analysis Baselines. We compare the proposed method with standard Cox models, f-DRO and its exact Cox variant (Exact-Cox) for latent subpopulation shift, and Robust-Cox (Minder and Bednarski, 1996) for outlier contamination. Both f-DRO and Exact-Cox use the χ2χ^2-divergence following Hu and Chen (2024); f-DRO uses sample splitting, whereas Exact-Cox uses the full likelihood without sample splitting. For nonlinear Robust-Cox, we extend the original estimator to a neural predictor. Datasets. We evaluate METABRIC (Curtis et al., 2012) and FLC (Dispenzieri et al., 2012), defining four subgroups in each dataset based on clinical covariates following Hu and Chen (2024). Group labels are unavailable during training and validation. Further details, including sample sizes, censoring rates, and subgroup definitions, are provided in Appendix A.3.5. Table 2: Results on FLC under subpopulation shift with outlier contamination. Settings and notation follow Table 1. Model Method IBS ↓ Test log-likelihood ↑ Overall C-index ↑ Worst-group C-index ↑ Test tail loss ↓ Linear (10%) Elastic-Cox 0.100 (0.003) -6.573 (0.066) 0.793 (0.003) 0.618 (0.038) 7.650 (0.114) Robust-Cox 0.096 (0.001) -6.491 (0.016) 0.793 (0.005) 0.611 (0.040) 7.854 (0.031) f-DRO 0.107 (0.003) -6.708 (0.031) 0.797 (0.004) 0.626 (0.037) 7.409 (0.035) Exact-Cox 0.103 (0.005) -6.634 (0.063) 0.789 (0.006) 0.609 (0.037) 7.581 (0.127) Ours 0.119 (0.002) -6.941 (0.009) 0.798 (0.004) 0.631 (0.036) 7.283 (0.006) Linear (5%) Elastic-Cox 0.100 (0.003) -6.573 (0.066) 0.793 (0.003) 0.618 (0.038) 7.650 (0.114) Robust-Cox 0.096 (0.001) -6.491 (0.016) 0.793 (0.005) 0.611 (0.040) 7.854 (0.031) f-DRO 0.125 (0.002) -7.052 (0.008) 0.795 (0.005) 0.626 (0.030) 7.294 (0.004) Exact-Cox 0.109 (0.010) -6.747 (0.192) 0.722 (0.140) 0.576 (0.085) 7.535 (0.072) Ours 0.119 (0.002) -6.941 (0.009) 0.798 (0.004) 0.631 (0.036) 7.283 (0.006) Nonlinear (10%) Cox-N 0.104 (0.006) -6.617 (0.065) 0.782 (0.010) 0.623 (0.029) 7.801 (0.208) Robust-Cox 0.106 (0.007) -6.690 (0.126) 0.760 (0.031) 0.602 (0.030) 8.134 (0.324) f-DRO 0.111 (0.004) -6.725 (0.033) 0.785 (0.009) 0.624 (0.028) 7.447 (0.043) Exact-Cox 0.112 (0.008) -6.783 (0.155) 0.723 (0.125) 0.561 (0.081) 7.622 (0.188) Ours 0.116 (0.003) -6.837 (0.057) 0.777 (0.019) 0.630 (0.033) 7.371 (0.067) Nonlinear (5%) Cox-N 0.104 (0.006) -6.617 (0.065) 0.782 (0.010) 0.623 (0.029) 7.801 (0.208) Robust-Cox 0.106 (0.007) -6.687 (0.129) 0.761 (0.032) 0.603 (0.031) 8.134 (0.324) f-DRO 0.126 (0.002) -7.055 (0.009) 0.779 (0.013) 0.620 (0.023) 7.298 (0.004) Exact-Cox 0.119 (0.012) -6.948 (0.214) 0.625 (0.181) 0.508 (0.117) 7.579 (0.207) Ours 0.123 (0.004) -6.976 (0.043) 0.781 (0.017) 0.628 (0.041) 7.284 (0.008) Table 3: Performance across hidden-layer sizes h on METABRIC under outlier contamination. Settings and notation follow the 10% validation setting in Table 1, whose nonlinear block corresponds to h=2h=2. h Method IBS ↓ Test log-likelihood ↑ Overall C-index ↑ Worst-group C-index ↑ Test tail loss ↓ 4 Cox-N 0.183 (0.007) -5.118 (0.047) 0.597 (0.012) 0.474 (0.073) 6.658 (0.188) Robust-Cox 0.183 (0.009) -5.139 (0.062) 0.612 (0.019) 0.506 (0.065) 6.914 (0.322) f-DRO 0.185 (0.004) -5.135 (0.021) 0.627 (0.028) 0.519 (0.036) 5.832 (0.021) Exact-Cox 0.199 (0.006) -5.225 (0.036) 0.452 (0.060) 0.391 (0.054) 6.245 (0.079) Ours 0.177 (0.004) -5.077 (0.023) 0.646 (0.027) 0.546 (0.042) 5.946 (0.045) 8 Cox-N 0.189 (0.005) -5.184 (0.068) 0.581 (0.029) 0.482 (0.043) 6.600 (0.662) Robust-Cox 0.186 (0.007) -5.179 (0.066) 0.611 (0.017) 0.504 (0.061) 7.106 (0.295) f-DRO 0.185 (0.005) -5.138 (0.037) 0.611 (0.012) 0.492 (0.055) 5.870 (0.053) Exact-Cox 0.201 (0.011) -5.243 (0.070) 0.439 (0.072) 0.387 (0.073) 6.344 (0.178) Ours 0.175 (0.005) -5.068 (0.037) 0.654 (0.027) 0.582 (0.045) 5.971 (0.084) 16 Cox-N 0.191 (0.006) -5.191 (0.070) 0.571 (0.032) 0.467 (0.074) 6.313 (0.731) Robust-Cox 0.190 (0.009) -5.231 (0.076) 0.608 (0.021) 0.477 (0.051) 7.320 (0.374) f-DRO 0.186 (0.004) -5.140 (0.028) 0.600 (0.020) 0.494 (0.056) 5.888 (0.049) Exact-Cox 0.201 (0.010) -5.250 (0.062) 0.465 (0.064) 0.394 (0.045) 6.367 (0.207) Ours 0.175 (0.003) -5.065 (0.028) 0.653 (0.017) 0.568 (0.057) 5.980 (0.089) Experimental Setup. Although the framework applies to general survival models, we consider linear Cox models and neural Cox models (Cox-N) with a single hidden layer of h=2h=2 units. Motivated by Duchi and Namkoong (2021), we select hyperparameters for Robust-Cox and the DRO-based methods by minimizing validation tail loss without using group labels. We evaluate overall performance using the integrated Brier score (IBS), test log-likelihood, and C-index, and robustness using worst-group C-index and 5% test tail loss. Details of the datasets, data splits, baselines, training, and hyperparameter selection are provided in Appendices A.3.5–A.3.7. Outlier Construction. We evaluate subpopulation shift with outlier contamination; results without outliers appear in Appendix A.3.8. We generate new survival times by inverting the fitted Cox risks of the highest-risk training samples: the top 2% for linear models and 5% for nonlinear models. These settings are motivated by survival outliers whose event times are unusually long or short relative to predictions based on their covariates (Luo et al., 2022). We use risk inversion to introduce such departures in a small fraction of the training samples, with effects that may vary across datasets and model classes. Construction details and additional contamination experiments are provided in Appendices A.3.5, A.3.9, and A.3.10. Results under outlier contamination. On METABRIC (Table 1), the proposed method outperforms all baselines in both overall and worst-group C-index in every setting, while achieving the second-lowest test tail loss. On FLC (Table 3), the results show a trade-off between average predictive performance and robustness: the proposed method achieves the highest worst-group C-index and the lowest test tail loss in every setting, but has higher IBS and lower test log-likelihood than the standard Cox models and Robust-Cox. Additionally, we conduct experiments with appended outliers on FLC using linear models and with 2% and 5% contamination across all dataset–model combinations (Appendices A.3.9 and A.3.10). Across these experiments, the proposed method consistently achieves the lowest or second-lowest test tail loss. Runtime comparisons are provided in Appendix A.3.11. Ablation across network configurations. We evaluate larger networks on METABRIC using hidden-layer sizes h=4,8,16h=4,8,16 under the same contamination and model-selection settings. As shown in Table 3, the proposed method achieves the highest overall and worst-group C-index, the lowest IBS, and the highest test log-likelihood in all three configurations. Its test tail loss is consistently second lowest, with f-DRO attaining the lowest values. These results show that the gains in overall and worst-group discrimination persist across the larger networks evaluated here. 7 Conclusion We proposed a coupled DRO framework for survival analysis that jointly addresses subpopulation shift and outlier contamination. The formulation couples nominal refinement with adversarial reweighting while preserving cross-sample interactions in non-decomposable survival losses. We derived outer updates from the inner problem’s KKT multipliers, yielding an alternating algorithm that directly optimizes the coupled objective. Experiments demonstrate improved worst-group C-index and consistently low test tail loss under contamination, with competitive overall C-index. Future work could develop ambiguity sets that also capture within-group covariate shifts while retaining outlier robustness. AI use statement We used generative AI tools to assist with manuscript drafting and editing and implementation details. The authors reviewed all AI-assisted content and take responsibility for the final manuscript. References Andreux et al. (2020) M. Andreux, A. Manoel, R. Menuet, C. Saillard, and C. Simpson Federated survival analysis with discrete-time Cox models. arXiv:2006.08997. Cited by: §2. Archetti et al. (2023) A. Archetti, E. Lomurno, F. Lattari, A. Martin, and M. Matteucci Heterogeneous datasets for federated survival analysis simulation. In International Conference on Performance Engineering Companion, p. 173–180. Cited by: §2. Archetti and Matteucci (2023) A. Archetti and M. Matteucci Federated survival forests. In International Joint Conference on Neural Networks, p. 1–9. Cited by: §2. Avati et al. (2020) A. Avati, T. Duan, S. Zhou, K. Jung, N. H. Shah, and A. Y. Ng Countdown regression: sharp and calibrated survival predictions. In Conference on Uncertainty in Artificial Intelligence, p. 145–155. Cited by: §2. Beck and Teboulle (2003) A. Beck and M. Teboulle Mirror descent and nonlinear projected subgradient methods for convex optimization. Operations Research Letters 31 (3), p. 167–175. Cited by: §A.1.4. Bennett (1983) S. Bennett Analysis of survival data by the proportional odds model. Statistics in Medicine 2 (2), p. 273–277. Cited by: §2. Bertsekas (1976) D. P. Bertsekas On the Goldstein–Levitin–Polyak gradient projection method. IEEE Transactions on Automatic Control 21 (2), p. 174–184. Cited by: §A.1.4. Buckley and James (1979) J. Buckley and I. James Linear regression with censored data. Biometrika 66 (3), p. 429–436. Cited by: §2. Cox (1972) D. R. Cox Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methodological) 34 (2), p. 187–202. Cited by: §1, §2, §4.3. Curtis et al. (2012) C. Curtis, S. P. Shah, S. Chin, G. Turashvili, O. M. Rueda, M. J. Dunning, D. Speed, A. G. Lynch, S. Samarajiwa, Y. Yuan, et al. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature 486 (7403), p. 346–352. Cited by: §A.3.5, §6.2. Delage and Ye (2010) E. Delage and Y. Ye Distributionally robust optimization under moment uncertainty with application to data-driven problems. Operations Research 58 (3), p. 595–612. Cited by: §2. Dispenzieri et al. (2012) A. Dispenzieri, J. A. Katzmann, R. A. Kyle, D. R. Larson, T. M. Therneau, C. L. Colby, R. J. Clark, G. P. Mead, S. Kumar, L. J. Melton, et al. Use of nonclonal serum immunoglobulin free light chains to predict overall survival in the general population. Mayo Clinic Proceedings 87 (6), p. 517–523. Cited by: §A.3.5, §6.2. Do et al. (2023) H. Do, Y. Chang, Y. S. Cho, P. Smyth, and J. Zhong Fair survival time prediction via mutual information minimization. In Machine Learning for Healthcare Conference, p. 128–149. Cited by: §2. Duchi and Namkoong (2021) J. C. Duchi and H. Namkoong Learning models with uniform performance via distributionally robust optimization. The Annals of Statistics 49 (3), p. 1378–1406. Cited by: §A.3.7, §1, §2, §3, §6.2. Fan et al. (2024) S. Fan, R. Xu, Q. Dong, Y. He, C. Chang, and P. Cui Stable Cox regression for survival analysis under distribution shifts. Nature Machine Intelligence 6 (12), p. 1525–1541. Cited by: §1. Farcomeni and Ventura (2012) A. Farcomeni and L. Ventura An overview of robust methods in medical research. Statistical Methods in Medical Research 21 (2), p. 111–133. Cited by: §2. Farcomeni and Viviani (2011) A. Farcomeni and S. Viviani Robust estimation for the Cox regression model based on trimming. Biometrical Journal 53 (6), p. 956–973. Cited by: §2. Gao et al. (2024) R. Gao, X. Chen, and A. J. Kleywegt Wasserstein distributionally robust optimization and variation regularization. Operations Research 72 (3), p. 1177–1191. Cited by: §2. Hashimoto et al. (2018) T. B. Hashimoto, M. Srivastava, H. Namkoong, and P. Liang Fairness without demographics in repeated loss minimization. In International Conference on Machine Learning, p. 1929–1938. Cited by: §2. Heritier et al. (2009) S. Heritier, E. Cantoni, S. Copt, and M. Victoria-Feser Robust methods in biostatistics. John Wiley & Sons. Cited by: §2. Hong et al. (2024) F. Hong, J. Yao, Y. Lyu, Z. Zhou, I. W. Tsang, Y. Zhang, and Y. Wang On harmonizing implicit subpopulations. In International Conference on Learning Representations, Cited by: §2. Hu and Chen (2022) S. Hu and G. H. Chen Distributionally robust survival analysis: a novel fairness loss without demographics. In Machine Learning for Health Symposium, p. 62–87. Cited by: §1, §2. Hu and Chen (2024) S. Hu and G. H. Chen Fairness in survival analysis with distributionally robust optimization. Journal of Machine Learning Research 25 (246), p. 1–85. Cited by: 3rd item, 4th item, §A.3.5, §1, §1, §1, §2, §4.3, §6.1, §6.2, §6.2. Huang et al. (2026) J. Huang, J. Lei, and G. Shi Robust learning via nested distributionally robust optimization. In International Conference on Machine Learning, Cited by: §2. Huber (1992) P. J. Huber Robust estimation of a location parameter. In Breakthroughs in statistics: Methodology and distribution, p. 492–518. Cited by: §3. Ishwaran et al. (2008) H. Ishwaran, U. B. Kogalur, E. H. Blackstone, and M. S. Lauer Random survival forests. The Annals of Applied Statistics 2 (3), p. 841–860. Cited by: §2. Jin et al. (2025) Y. Jin, L. Wise, and I. Ch. Paschalidis Distributionally robust learning in survival analysis. In Conference on Health, Inference, and Learning, p. 369–380. Cited by: §1, §4.3. Jo et al. (2026) S. H. Jo, S. Kim, and M. Chae Mitigating spurious correlation via distributionally robust learning with hierarchical ambiguity sets. In International Conference on Learning Representations, Cited by: §2. Kalbfleisch and Prentice (2002) J. D. Kalbfleisch and R. L. Prentice The statistical analysis of failure time data. 2nd edition, John Wiley & Sons. Cited by: §2. Kaplan and Meier (1958) E. L. Kaplan and P. Meier Nonparametric estimation from incomplete observations. Journal of the American Statistical Association 53 (282), p. 457–481. Cited by: §2. Katzman et al. (2018) J. L. Katzman, U. Shaham, A. Cloninger, J. Bates, T. Jiang, and Y. Kluger DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network. BMC Medical Research Methodology 18, p. 24. Cited by: §2. Keya et al. (2021) K. N. Keya, R. Islam, S. Pan, I. Stockwell, and J. Foulds Equitable allocation of healthcare resources with fair survival models. In SIAM International Conference on Data Mining, p. 190–198. Cited by: §2. Kim et al. (2026) S. Kim, S. H. Jo, W. Ha, and M. Chae Distributionally robust classification for multi-source unsupervised domain adaptation. In International Conference on Learning Representations, Cited by: §2. Kuhn et al. (2019) D. Kuhn, P. Mohajerin Esfahani, V. A. Nguyen, and S. Shafieezadeh-Abadeh Wasserstein distributionally robust optimization: theory and applications in machine learning. In Operations Research & Management Science in the Age of Analytics, p. 130–166. Cited by: §2. Kvamme et al. (2019) H. Kvamme, Ø. Borgan, and I. Scheel Time-to-event prediction with neural networks and Cox regression. Journal of Machine Learning Research 20 (129), p. 1–30. Cited by: §5. Lee et al. (2018) C. Lee, W. R. Zame, J. Yoon, and M. van der Schaar DeepHit: a deep learning approach to survival analysis with competing risks. In AAAI Conference on Artificial Intelligence, p. 2314–2321. Cited by: §1, §2. Lin et al. (2022) F. Lin, X. Fang, and Z. Gao Distributionally robust optimization: a review on theory and applications. Numerical Algebra, Control and Optimization 12 (1), p. 159–212. Cited by: §2. Liu et al. (2021) E. Z. Liu, B. Haghgoo, A. S. Chen, A. Raghunathan, P. W. Koh, S. Sagawa, P. Liang, and C. Finn Just train twice: improving group robustness without training group information. In International Conference on Machine Learning, p. 6781–6792. Cited by: §2. Liu et al. (2014) Z. Liu, X. Zhang, and S. Zhang Breast tumor subgroups reveal diverse clinical prognostic power. Scientific Reports 4, p. 4002. Cited by: §A.3.5. Luo et al. (2022) B. Luo, X. Gao, and S. Halabi Penalized weighted proportional hazards model for robust variable selection and outlier detection. Statistics in Medicine 41 (17), p. 3398–3420. Cited by: §A.3.8, §2, §6.2. Minder and Bednarski (1996) C. E. Minder and T. Bednarski A robust method for proportional hazards regression. Statistics in Medicine 15 (10), p. 1033–1047. Cited by: 5th item, §6.2. Mohajerin Esfahani and Kuhn (2018) P. Mohajerin Esfahani and D. Kuhn Data-driven distributionally robust optimization using the Wasserstein metric: performance guarantees and tractable reformulations. Mathematical Programming 171 (1–2), p. 115–166. Cited by: §1, §2. Namkoong and Duchi (2016) H. Namkoong and J. C. Duchi Stochastic gradient methods for distributionally robust optimization with f-divergences. In Advances in Neural Information Processing Systems, Cited by: §2, §3. Rahman and Purushotham (2022) M. M. Rahman and S. Purushotham Fair and interpretable models for survival analysis. In SIGKDD Conference on Knowledge Discovery and Data Mining, p. 1452–1462. Cited by: §1, §2. Rahman and Purushotham (2023) M. M. Rahman and S. Purushotham FedPseudo: privacy-preserving pseudo value-based deep learning models for federated survival analysis. In SIGKDD Conference on Knowledge Discovery and Data Mining, p. 1999–2009. Cited by: §2. Rahman and Purushotham (2025) M. M. Rahman and S. Purushotham Fair federated survival analysis. In AAAI Conference on Artificial Intelligence, p. 20104–20112. Cited by: §2. Ranganath et al. (2016) R. Ranganath, A. Perotte, N. Elhadad, and D. Blei Deep survival analysis. In Machine Learning for Healthcare Conference, p. 101–114. Cited by: §2. Raykar et al. (2007) V. C. Raykar, H. Steck, B. Krishnapuram, C. Dehing-Oberije, and P. Lambin On ranking in survival analysis: bounds on the concordance index. In Advances in Neural Information Processing Systems, p. 1209–1216. Cited by: §1. Rockafellar and Wets (1998) R. T. Rockafellar and R. J. B. Wets Variational analysis. Springer. Cited by: §A.1.2. Rudner et al. (2024) T. G. J. Rudner, Y. S. Zhang, A. G. Wilson, and J. Kempe Mind the GAP: improving robustness to subpopulation shifts with group-aware priors. In International Conference on Artificial Intelligence and Statistics, p. 127–135. Cited by: §2. Sagawa et al. (2020) S. Sagawa, P. W. Koh, T. B. Hashimoto, and P. Liang Distributionally robust neural networks for group shifts: on the importance of regularization for worst-case generalization. In International Conference on Learning Representations, Cited by: §1, §2, §3, §5, §6.1. Seidi et al. (2024) N. Seidi, S. Roy, S. K. Das, and A. Tripathy Addressing data heterogeneity in federated learning of Cox proportional hazards models. In IEEE International Conference on E-health Networking, Application & Services, Cited by: §2. Shapiro et al. (2021) A. Shapiro, D. Dentcheva, and A. Ruszczynski Lectures on stochastic programming: Modeling and theory. 3rd edition, SIAM. Cited by: §A.1.1, §2, §4.1. Sinha (2019) S. K. Sinha Robust estimation in accelerated failure time models. Lifetime Data Analysis 25 (1), p. 52–78. Cited by: §2. Stromberg et al. (2024) N. Stromberg, R. Ayyagari, M. Welfert, S. Koyejo, R. Nock, and L. Sankar Robustness to subpopulation shift with domain label noise via regularized annotation of domains. Transactions on Machine Learning Research. Cited by: §2. Tang et al. (2023) W. Tang, K. He, G. Xu, and J. Zhu Survival analysis via ordinary differential equations. Journal of the American Statistical Association 118 (544), p. 2406–2421. Cited by: §2. Wang et al. (2022) X. Wang, H. G. Zhang, X. Xiong, C. Hong, G. M. Weber, G. A. Brat, C. Bonzel, Y. Luo, R. Duan, N. P. Palmer, et al. SurvMaximin: Robust federated approach to transporting survival risk prediction models. Journal of Biomedical Informatics 134, p. 104176. Cited by: §2. Yang et al. (2023) Y. Yang, H. Zhang, D. Katabi, and M. Ghassemi Change is hard: a closer look at subpopulation shift. In International Conference on Machine Learning, p. 39584–39622. Cited by: §1. Zeng et al. (2026) L. Zeng, W. Tang, Z. Ren, and Y. Ding Mini-batch estimation for deep Cox models: Statistical foundations and practical guidance. Journal of the American Statistical Association 121 (554), p. 988–999. Cited by: §5. Zeng et al. (2023) Y. Zeng, K. Greenewald, L. Jung, K. Lee, J. Solomon, and M. Yurochkin Outlier-robust group inference via gradient space clustering. In NeurIPS Workshop on Distribution Shifts, Cited by: §2. Zhai et al. (2021) R. Zhai, C. Dan, Z. Kolter, and P. Ravikumar DORO: Distributional and outlier robust optimization. In International Conference on Machine Learning, p. 12345–12355. Cited by: §1, §2, §3, §4.2. Zhang and Weiss (2022) W. Zhang and J. C. Weiss Longitudinal fairness with censorship. In AAAI Conference on Artificial Intelligence, p. 12235–12243. Cited by: §2. Zhang and Weiss (2023) W. Zhang and J. C. Weiss Fairness with censorship and group constraints. Knowledge and Information Systems 65 (6), p. 2571–2594. Cited by: §2. Appendix A Appendix Contents of the appendix A.1 Formulation and algorithm details For additive losses, ℒ(q,θ)=∑i=1nqiℓ(zi,θ)L(q,θ)= _i=1^nq_i (z_i,θ) is linear in q, so the optimal weights follow the ordering of the losses. We first derive this trimming and reweighting structure, then use KKT multipliers to obtain gradient-based updates for non-decomposable losses such as the weighted Cox loss. We also analyze computational complexity and local convergence. For fixed θ, recall the inner value function Φα(q′,θ):=maxq∈α(q′)ℒ(q,θ). _α(q ;θ):= _q _α(q )L(q,θ). (9) We omit θ when it is fixed. A.1.1 Additive losses Write ℓi:=ℓ(zi,θ) _i:= (z_i,θ) and Φα(q′):=Φα(q′,θ) _α(q ):= _α(q ;θ). We first express the ambiguity set as bounds on sample weights, then derive ϵε-trimming for α=1α=1 and the two-stage solution for general α. The uniform-weight case gives the explicit allocation in Corollary 1. Lemma 1 (Equivalent discrete representation of the ambiguity set). Let P¯ P be the nominal distribution supported on D, with corresponding weight vector p¯∈n−1 p _n-1. For γ∈(0,1]γ∈(0,1], the set γ(P¯)Q_γ( P) admits the equivalent finite-support representation (5), i.e., γ(p¯)=q∈n−1:0≤qi≤p¯iγ,∀i∈[n].Q_γ( p)= \q _n-1:0≤ q_i≤ p_iγ,\;∀ i∈[n] \. Moreover, every such q can be reparameterized as qi=p¯iγwi,w∈[0,1]n,∑i=1np¯iwi=γ. q_i= p_iγw_i, w∈[0,1]^n, _i=1^n p_iw_i=γ. (10) Proof. We first consider the case γ∈(0,1)γ∈(0,1). (⇒ ) Suppose that Q∈γ(P¯)Q _γ( P). Since P¯ P is supported on the dataset =zii=1nD=\z_i\_i=1^n, there exists a weight vector p¯∈n−1 p _n-1 such that P¯=∑i=1np¯iδzi, P= _i=1^n p_i _z_i, where δzi _z_i denotes the Dirac measure at ziz_i. By the definition of γ(P¯)Q_γ( P), there exists a probability distribution Q~∈() Q (Z) such that P¯=γQ+(1−γ)Q~. P=γ Q+(1-γ) Q. Since P¯ P is supported on D, the above identity implies that Q and Q~ Q are also supported on D. Hence there exist weight vectors q,q~∈n−1q, q _n-1 such that Q=∑i=1nqiδzi,Q~=∑i=1nq~iδzi.Q= _i=1^nq_i _z_i, Q= _i=1^n q_i _z_i. By comparing the coefficients of the Dirac measures at each ziz_i, we obtain p¯=γq+(1−γ)q~. p=γ q+(1-γ) q. Taking the i-th coordinate and using q~i≥0 q_i≥ 0, we have p¯i=γqi+(1−γ)q~i≥γqi, p_i=γ q_i+(1-γ) q_i≥γ q_i, which implies qi≤p¯iγ.q_i≤ p_iγ. Since q∈n−1q _n-1, we also have qi≥0q_i≥ 0 for all i and ∑i=1nqi=1 _i=1^nq_i=1. Hence q∈γ(p¯)q _γ( p) in (5). (⇐ ) Conversely, suppose that q∈n−1q _n-1 satisfies 0≤qi≤p¯iγ,∀i∈[n].0≤ q_i≤ p_iγ, ∀ i∈[n]. Define a vector q~∈ℝn q ^n by q~i:=p¯i−γqi1−γ,i∈[n]. q_i:= p_i-γ q_i1-γ, i∈[n]. Since qi≤p¯i/γq_i≤ p_i/γ, we have q~i≥0 q_i≥ 0 for all i∈[n]i∈[n]. Moreover, ∑i=1nq~i=11−γ∑i=1n(p¯i−γqi)=11−γ(1−γ⋅1)=1. _i=1^n q_i= 11-γ _i=1^n ( p_i-γ q_i )= 11-γ (1-γ· 1 )=1. Hence q~∈n−1 q _n-1. Now define probability distributions Q=∑i=1nqiδzi,Q~=∑i=1nq~iδzi.Q= _i=1^nq_i _z_i, Q= _i=1^n q_i _z_i. Then γQ+(1−γ)Q~=∑i=1n(γqi+(1−γ)q~i)δzi=∑i=1np¯iδzi=P¯.γ Q+(1-γ) Q= _i=1^n (γ q_i+(1-γ) q_i ) _z_i= _i=1^n p_i _z_i= P. Therefore, Q∈γ(P¯)Q _γ( P). We now show the reparameterization in (10). Suppose that q satisfies (5). For each i∈[n]i∈[n] with p¯i>0 p_i>0, define wi:=γqip¯i.w_i:= γ q_i p_i. If p¯i=0 p_i=0, then (5) implies qi=0q_i=0, and we define wi:=0w_i:=0. Then wi≥0w_i≥ 0 for all i. Moreover, from qi≤p¯i/γq_i≤ p_i/γ, it follows that wi≤1w_i≤ 1 whenever p¯i>0 p_i>0. Hence w∈[0,1]nw∈[0,1]^n. By construction, qi=p¯iγwi,∀i∈[n].q_i= p_iγw_i, ∀ i∈[n]. Finally, ∑i=1np¯iwi=∑i:p¯i>0p¯iγqip¯i=γ∑i=1nqi=γ, _i=1^n p_iw_i= _i: p_i>0 p_i γ q_i p_i=γ _i=1^nq_i=γ, since q∈n−1q _n-1. If γ=1γ=1, then the condition P¯=Q P=Q implies q=p¯q= p. The discrete constraint becomes 0≤qi≤p¯i0≤ q_i≤ p_i with q∈n−1q _n-1, which also implies q=p¯q= p. The reparameterization holds with wi=1w_i=1 for p¯i>0 p_i>0 and wi=0w_i=0 for p¯i=0 p_i=0. This completes the proof. ∎ Using the representation in Lemma 1, we first consider the special case α=1α=1. In this case, the inner ambiguity set α(q′)Q_α(q ) reduces to the singleton q′\q \, and the coupled DRO problem reduces to min∑i=1nq′∈1−ϵ(p¯)qi′ℓi. _q _1-ε( p) _i=1^nq _i _i. (11) The minimizer assigns mass to smaller-loss samples first, yielding ϵε-trimming. Proposition 2 (Exact ϵε-trimming under α=1α=1). Assume α=1α=1. Consider the problem q′⋆∈argminq′∈1−ϵ(p¯)∑i=1nq′iℓi. q ∈ _q _1-ε( p) _i=1^nq _i _i. Then there exists an optimal solution of the form qi′⋆=p¯iwi⋆1−ϵ, q_i = p_iw_i 1-ε, where w⋆∈[0,1]nw ∈[0,1]^n satisfies ∑i=1np¯iwi⋆=1−ϵ, _i=1^n p_iw_i =1-ε, and, after reindexing so that ℓ1≤ℓ2≤⋯≤ℓn, _1≤ _2≤·s≤ _n, there exists k∈[n]k∈[n] such that wi⋆=1,i<k,γ,i=k,0,i>k,w_i = cases1,&i<k,\\ γ,&i=k,\\ 0,&i>k, cases for some γ∈[0,1]γ∈[0,1] chosen to satisfy the probability mass constraint ∑i=1nqi′∗=1 _i=1^nq *_i=1. Proof. The argument transfers mass from larger to smaller losses without increasing the objective. By Lemma 1, any q′∈1−ϵ(p¯)q _1-ε( p) admits the representation qi′=p¯iwi1−ϵ,w∈[0,1]n,∑i=1np¯iwi=1−ϵ.q _i= p_iw_i1-ε, w∈[0,1]^n, _i=1^n p_iw_i=1-ε. Then, (11) reduces to minw∈[0,1]n _w∈[0,1]^n ∑i=1np¯iwiℓis.t.∑i=1np¯iwi=1−ϵ. _i=1^n p_iw_i _i .t. _i=1^n p_iw_i=1-ε. Let mi:=p¯iwim_i:= p_iw_i. The problem is equivalently min∑i=1n0≤mi≤p¯iℓimis.t.∑i=1nmi=1−ϵ. _0≤ m_i≤ p_i _i=1^n _im_i .t. _i=1^nm_i=1-ε. The objective is linear in the allocated mass mim_i with coefficient ℓi _i, while p¯i p_i acts as an upper bound. Reindex so that ℓ1≤ℓ2≤⋯≤ℓn _1≤ _2≤·s≤ _n. Suppose there exist i<ji<j such that mi<p¯im_i< p_i and mj>0m_j>0. Then for any 0<δ≤minp¯i−mi,mj0<δ≤ \ p_i-m_i,\;m_j\, the perturbation mi←mi+δ,mj←mj−δm_i← m_i+δ, m_j← m_j-δ preserves feasibility and changes the objective by δ(ℓi−ℓj)≤0δ( _i- _j)≤ 0, with strict inequality if ℓi<ℓj _i< _j. Hence an optimal solution assigns mass in increasing loss order up to the bounds, including when losses are tied. Thus there exists an index k∈[n]k∈[n] such that mi⋆=p¯i,i<k,m~k,i=k,0,i>k,m_i = cases p_i,&i<k,\\ m_k,&i=k,\\ 0,&i>k, cases for some m~k∈[0,p¯k] m_k∈[0, p_k] chosen to satisfy the probability mass constraint. Equivalently, there exist k∈[n]k∈[n] and γ∈[0,1]γ∈[0,1] such that wi⋆=1,i<k,γ,i=k,0,i>k,w_i = cases1,&i<k,\\ γ,&i=k,\\ 0,&i>k, cases with ∑i=1np¯iwi⋆=1−ϵ _i=1^n p_iw_i =1-ε. This completes the proof. ∎ The next result extends this structure to the general case α∈(0,1]α∈(0,1] under ϵ∈[0,1)ε∈[0,1). We first show that the inner maximization problem Φα(q′) _α(q ) admits a weighted CVaR representation. Lemma 2 (Weighted CVaR dual representation of Φα _α). For any q′∈n−1q _n-1, the inner maximization problem Φα(q′) _α(q ) admits the dual representation Φα(q′)=minτ∈ℝτ+1α∑i=1nqi′(ℓi−τ)+, _α(q )= _τ \τ+ 1α _i=1^nq_i ( _i-τ)_+ \, where (x)+:=maxx,0(x)_+:= \x,0\. Proof. Recall that the inner maximization problem is Φα(q′)=max∑i=1nqiℓi:q∈n−1, 0≤qi≤1αqi′,∀i∈[n]. _α(q )= \ _i=1^nq_i _i:q _n-1,\ 0≤ q_i≤ 1αq_i ,\ ∀ i∈[n] \. This is precisely the standard primal formulation of weighted CVaR under the reference weight vector q′q . By the dual representation of CVaR (Shapiro et al., 2021), we obtain Φα(q′)=minτ∈ℝτ+1α∑i=1nqi′(ℓi−τ)+. _α(q )= _τ \τ+ 1α _i=1^nq_i ( _i-τ)_+ \. This completes the proof. ∎ The weighted CVaR representation yields the following two-stage solution. Theorem 1 (Two-stage structure of the coupled DRO problem). There exists an optimal solution with the following structure. After reordering the indices so that ℓ1≤ℓ2≤⋯≤ℓn, _1≤ _2≤·s≤ _n, i.e., in increasing order of losses, the solution exhibits a two-stage structure consisting of an outer and an inner component. (Outer structure) The outer minimization assigns mass to smaller-loss samples first. Specifically, qi′⋆=p¯iwi⋆1−ϵ,w⋆∈[0,1]n,∑i=1np¯iwi⋆=1−ϵ,q_i = p_iw_i 1-ε, w ∈[0,1]^n, _i=1^n p_iw_i =1-ε, where wi⋆=1,i<k,γ,i=k,0,i>k,w_i = cases1,&i<k,\\ γ,&i=k,\\ 0,&i>k, cases for some k∈[n]k∈[n] and γ∈[0,1]γ∈[0,1]. (Inner structure) Given q′⋆q , the inner maximization concentrates mass on larger-loss samples. In particular, qi⋆=0,i<t,η1αqt′⋆,i=t,1αqi′⋆,i>t,∑i=1nqi⋆=1,q_i = cases0,&i<t,\\ η\, 1αq_t ,&i=t,\\ 1αq_i ,&i>t, cases _i=1^nq_i =1, for some t∈[n]t∈[n] and η∈[0,1]η∈[0,1]. Proof. The CVaR representation turns the outer problem into a joint minimization over q′q and a threshold. Mass transfers then characterize both the outer and inner optimizers. By Lemma 2, for any q′∈n−1q _n-1, Φα(q′)=minτ∈ℝτ+1α∑i=1nqi′(ℓi−τ)+. _α(q )= _τ \τ+ 1α _i=1^nq _i( _i-τ)_+ \. Therefore, the DRO objective can be written as the joint minimization minq′∈1−ϵ(p¯)Φα(q′)=minq′∈1−ϵ(p¯),τ∈ℝτ+1α∑i=1nqi′(ℓi−τ)+. _q _1-ε( p) _α(q )= _q _1-ε( p),\ τ \τ+ 1α _i=1^nq _i( _i-τ)_+ \. (12) Let q′⋆q be an outer minimizer of Φα(q′) _α(q ), and choose τ⋆∈argminτ∈ℝτ+1α∑iqi′⋆(ℓi−τ)+.τ ∈ _τ \τ+ 1α _iq_i ( _i-τ)_+ \. For any q′∈1−ϵ(p¯)q _1-ε( p) and τ∈ℝτ , τ+1α∑i=1nqi′(ℓi−τ)+≥Φα(q′)≥Φα(q′⋆)=τ⋆+1α∑i=1nqi′⋆(ℓi−τ⋆)+.τ+ 1α _i=1^nq _i( _i-τ)_+≥ _α(q )≥ _α(q )=τ + 1α _i=1^nq_i ( _i-τ )_+. Hence (q′⋆,τ⋆)(q ,τ ) attains the minimum in (12). (Outer structure). Define ai⋆:=(ℓi−τ⋆)+a_i :=( _i-τ )_+. For fixed τ⋆τ , the objective is linear in q′q , so q′⋆q solves min∑i=1nq′∈1−ϵ(p¯)qi′ai⋆. _q _1-ε( p) _i=1^nq _ia_i . By Lemma 1, write qi′=p¯iwi/(1−ϵ)q _i= p_iw_i/(1-ε) with w∈[0,1]nw∈[0,1]^n and ∑ip¯iwi=1−ϵ _i p_iw_i=1-ε. Equivalently, letting mi:=p¯iwim_i:= p_iw_i, we obtain min∑i=1n0≤mi≤p¯iai⋆mis.t.∑i=1nmi=1−ϵ. _0≤ m_i≤ p_i _i=1^na_i m_i .t. _i=1^nm_i=1-ε. Reindex the samples so that ℓ1≤⋯≤ℓn _1≤·s≤ _n. Since ai⋆=(ℓi−τ⋆)+a_i =( _i-τ )_+ is non-decreasing in ℓi _i, this indexing also satisfies a1⋆≤⋯≤an⋆a_1 ≤·s≤ a_n . If there exist i<ji<j with mi<p¯im_i< p_i and mj>0m_j>0, then for any 0<δ≤minp¯i−mi,mj0<δ≤ \ p_i-m_i,\ m_j\, the perturbation mi←mi+δm_i← m_i+δ, mj←mj−δm_j← m_j-δ preserves feasibility and changes the objective by δ(ai⋆−aj⋆)≤0δ(a_i -a_j )≤ 0 with strict decrease if ai⋆<aj⋆a_i <a_j . Therefore, an optimal solution can be chosen that fills indices in this order, including when the ai⋆a_i are tied. Consequently, there exist k∈[n]k∈[n] and γ∈[0,1]γ∈[0,1] such that wi⋆=1,i<k,γ,i=k,0,i>k,w_i = cases1,&i<k,\\ γ,&i=k,\\ 0,&i>k, cases which shows that w⋆w has a threshold structure with at most one fractional component. (Inner structure). Given q′⋆q , the inner problem is the linear program max∑i=1nq∈α(q′⋆)qiℓi. _q _α(q ) _i=1^nq_i _i. Using the same loss ordering ℓ1≤⋯≤ℓn _1≤·s≤ _n, suppose that there exist i<ji<j such that qi>0q_i>0 and qj<1αqj′⋆q_j< 1αq_j . Then for any 0<δ≤minqi,1αqj′⋆−qj0<δ≤ \q_i,\ 1αq_j -q_j\, the perturbation qi←qi−δq_i← q_i-δ, qj←qj+δq_j← q_j+δ preserves feasibility and increases the objective by δ(ℓj−ℓi)≥0δ( _j- _i)≥ 0, with strict increase if ℓj>ℓi _j> _i. Therefore, an optimal solution can be chosen that fills indices in decreasing loss order, including when losses are tied, up to the caps qi≤1αqi′⋆q_i≤ 1αq_i . Consequently, there exist t∈[n]t∈[n] and η∈[0,1]η∈[0,1] such that qi⋆=0,i<t,η1αqt′⋆,i=t,1αqi′⋆,i>t,q_i = cases0,&i<t,\\ η\, 1αq_t ,&i=t,\\ 1αq_i ,&i>t, cases which shows that q⋆q has a threshold structure with at most one fractional boundary index. ∎ For uniform nominal weights, the two-stage solution takes the following explicit form. Corollary 1 (Exact trimming followed by exact worst-case allocation). Assume p¯i=1/n p_i=1/n for all i, and reorder the samples so that ℓ1≤ℓ2≤⋯≤ℓn. _1≤ _2≤·s≤ _n. Assume that (1−ϵ)n∈ℤ(1-ε)n and α(1−ϵ)n∈ℤα(1-ε)n . Let kϵ:=(1−ϵ)n,mα:=α(1−ϵ)n.k_ε:=(1-ε)n, m_α:=α(1-ε)n. Then there exists an optimal solution (q′⋆,q⋆)(q ,q ) such that qi′⋆=1(1−ϵ)n,i=1,…,kϵ,0,i>kϵ,q_i = cases 1(1-ε)n,&i=1,…,k_ε,\\[6.0pt] 0,&i>k_ε, cases and qi⋆=0,i=1,…,kϵ−mα,1α(1−ϵ)n,i=kϵ−mα+1,…,kϵ,0,i>kϵ.q_i = cases0,&i=1,…,k_ε-m_α,\\[6.0pt] 1α(1-ε)n,&i=k_ε-m_α+1,…,k_ε,\\[6.0pt] 0,&i>k_ε. cases In particular, q′⋆q is the exact ϵε-trimmed empirical distribution, and q⋆q places uniform mass on the largest α(1−ϵ)nα(1-ε)n losses among the retained samples. Proof. By Theorem 1, there exists an optimal solution (q′⋆,q⋆)(q ,q ) such that the outer optimizer is obtained by filling the smallest-loss indices first. Then, we have qi′⋆=p¯iwi⋆1−ϵ,q_i = p_iw_i 1-ε, where for some k and γ∈[0,1]γ∈[0,1], wi⋆=1,i<k,γ,i=k,0,i>k,∑i=1np¯iwi⋆=1−ϵ.w_i = cases1,&i<k,\\ γ,&i=k,\\ 0,&i>k, cases _i=1^n p_iw_i =1-ε. Since p¯i=1/n p_i=1/n, the normalization condition becomes ∑i=1np¯iwi⋆=1n((k−1)+γ)=1−ϵ. _i=1^n p_iw_i = 1n ((k-1)+γ )=1-ε. Because (1−ϵ)n=kϵ∈ℤ(1-ε)n=k_ε , it follows that (k−1)+γ=kϵ.(k-1)+γ=k_ε. Since kϵ∈ℤk_ε and kϵ=(k−1)+γk_ε=(k-1)+γ with γ∈[0,1]γ∈[0,1], we may take γ=1γ=1 and k=kϵk=k_ε, which also covers ϵ=0ε=0. Then, we have qi′⋆=1(1−ϵ)n,i=1,…,kϵ,0,i>kϵ.q_i = cases 1(1-ε)n,&i=1,…,k_ε,\\[6.0pt] 0,&i>k_ε. cases Now consider the inner maximization max∑i=1nq∈α(q′⋆)qiℓi,0≤qi≤1αqi′⋆,∑i=1nqi=1. _q _α(q ) _i=1^nq_i _i, 0≤ q_i≤ 1αq_i , _i=1^nq_i=1. Since ℓ1≤⋯≤ℓn _1≤·s≤ _n, the objective is maximized by assigning as much mass as possible to the largest-loss indices. Because the constraint qi≤1αqi′⋆q_i≤ 1αq_i holds for every i, and qi′⋆=0q_i =0 for all i>kϵi>k_ε, we have qi⋆=0for all i>kϵ.q_i =0 all i>k_ε. For each retained index i≤kϵi≤ k_ε, we have 1αqi′⋆=1α(1−ϵ)n=1mα, 1αq_i = 1α(1-ε)n= 1m_α, since mα=α(1−ϵ)nm_α=α(1-ε)n. Thus every non-zero upper bound is the same, namely 1/mα1/m_α. Therefore the maximizer is obtained by saturating exactly the largest mαm_α indices among 1,…,kϵ\1,…,k_ε\, namely i=kϵ−mα+1,…,kϵi=k_ε-m_α+1,…,k_ε. This yields qi⋆=0,i=1,…,kϵ−mα,1α(1−ϵ)n,i=kϵ−mα+1,…,kϵ,0,i>kϵ.q_i = cases0,&i=1,…,k_ε-m_α,\\[6.0pt] 1α(1-ε)n,&i=k_ε-m_α+1,…,k_ε,\\[6.0pt] 0,&i>k_ε. cases Since there are exactly mαm_α such indices, the total mass is mα⋅1α(1−ϵ)n=mαmα=1,m_α· 1α(1-ε)n= m_αm_α=1, so the simplex constraint is satisfied. This proves the claim. ∎ Remark. When (1−ϵ)n∉ℤ(1-ε)n , the loss-ordered outer solution has a fractional boundary atom at index ⌊(1−ϵ)n⌋+1 (1-ε)n +1. The inner optimizer fills the retained atoms in decreasing loss order, with partial allocation at the cutoff when needed. For uniform nominal weights, Algorithm 2 computes the weights by sorting the losses, trimming the largest losses, and emphasizing the largest retained losses. It assumes that (1−ϵ)n(1-ε)n and α(1−ϵ)nα(1-ε)n are integers; otherwise, fractional allocation is required at the boundaries as described above. Algorithm 2 Procedure for minimizing the coupled DRO objective under additive losses 1: data zii=1n\z_i\_i=1^n, parameters ϵ,αε,α, step size η 2: Initialize θ 3: Set kϵ←⌊(1−ϵ)n⌋k_ε← (1-ε)n and mα←⌊α(1−ϵ)n⌋m_α← α(1-ε)n 4: while not converged do 5: Compute losses ℓi=ℓ(zi,θ) _i= (z_i,θ) for i=1,…,ni=1,…,n 6: Reorder the indices so that ℓ1≤ℓ2≤⋯≤ℓn _1≤ _2≤·s≤ _n 7: Set qi′←1(1−ϵ)n,i=1,…,kϵ,0,i>kϵq_i ← cases 1(1-ε)n,&i=1,…,k_ε,\\ 0,&i>k_ε cases 8: Set qi←0,i=1,…,kϵ−mα,1α(1−ϵ)n,i=kϵ−mα+1,…,kϵ,0,i>kϵq_i← cases0,&i=1,…,k_ε-m_α,\\ 1α(1-ε)n,&i=k_ε-m_α+1,…,k_ε,\\ 0,&i>k_ε cases 9: Update θ←θ−η∇θ(∑i=1nqiℓ(zi,θ))θ←θ-η _θ ( _i=1^nq_i\, (z_i,θ) ) 10: end while 11: return θ A.1.2 Non-decomposable losses For non-decomposable losses, interactions across samples generally prevent the explicit weight allocation derived above. We instead alternate between ascent in q, an update of q′q based on an estimated supergradient, and descent in θ. The main step is to characterize the sensitivity of the inner value function Φα(q′,θ)=maxq∈α(q′)ℒ(q,θ). _α(q ;θ)= _q _α(q )L(q,θ). (13) Proposition 1 relates this sensitivity to the upper-bound KKT multipliers under local stability (Rockafellar and Wets, 1998). Theorem 2 gives their explicit form. Proposition 1 (restated). Fix θ, 0<α<10<α<1, and q′∈n−1q _n-1 with qi′>0q _i>0 for all i. Set b=q′/αb=q /α, ℱ(b)=q∈n−1:q≤bF(b)=\q _n-1:q≤ b\, and V(b):=maxq∈ℱ(b)ℒ(q,θ),Φα(q′,θ)=V(b).V(b):= _q (b)L(q,θ), _α(q ;θ)=V(b). Let q⋆(b)q (b) be a global maximizer, and assume that ℒ(⋅,θ)L(·,θ) is continuous on the simplex and continuously differentiable near q⋆(b)q (b). Assume local stability: for every sufficiently small Δb b, a global maximizer q⋆(b+Δb)q (b+ b) can be selected such that ‖q⋆(b+Δb)−q⋆(b)‖≤C‖Δb‖\|q (b+ b)-q (b)\|≤ C\| b\| (14) for some constant C>0C>0 independent of Δb b. Then any KKT multiplier μ⋆μ for the upper-bound constraints qi≤biq_i≤ b_i at q⋆(b)q (b) satisfies μ⋆α∈∂+Φα(q′,θ), μ α∈∂^+ _α(q ;θ), where ∂+∂^+ denotes the superdifferential: s∈∂+F(u)s∈∂^+F(u) if F(u+h)≤F(u)+⟨s,h⟩+o(‖h‖)F(u+h)≤ F(u)+ s,h +o(\|h\|) as h→0h→ 0 within the domain of F. Proof. We prove the result by combining the KKT conditions with a first-order Taylor expansion. The local stability assumption controls the remainder under perturbations of b. Write Δq=q⋆(b+Δb)−q⋆(b) q=q (b+ b)-q (b) and g⋆=∇qℒ(q⋆(b),θ)g = _qL(q (b),θ). The KKT conditions give g⋆=ν⋆−λ⋆+μ⋆g =ν 1-λ +μ , where λ⋆λ and μ⋆μ are the nonnegative lower- and upper-bound multipliers. Since both solutions sum to one, feasibility and complementary slackness imply ⟨g⋆,Δq⟩ g , q =−⟨λ⋆,Δq⟩+⟨μ⋆,Δq⟩ =- λ , q + μ , q ≤⟨μ⋆,Δb⟩. ≤ μ , b . By (14), ‖Δq‖=O(‖Δb‖)\| q\|=O(\| b\|). A first-order Taylor expansion therefore yields V(b+Δb)−V(b) V(b+ b)-V(b) =⟨g⋆,Δq⟩+o(‖Δb‖) = g , q +o(\| b\|) ≤⟨μ⋆,Δb⟩+o(‖Δb‖). ≤ μ , b +o(\| b\|). Thus μ⋆∈∂+V(b)μ ∈∂^+V(b). Substituting b=q′/αb=q /α gives μ⋆/α∈∂+Φα(q′,θ)μ /α∈∂^+ _α(q ;θ). ∎ Applying the supergradient result to weighted Cox. The following lemma gives conditions under which the weighted Cox loss satisfies Proposition 1. The key step is to bound changes in the inner optimizer by changes in its upper bounds. For fixed θ with fθ(xj)∈ℝf_θ(x_j) for all j, write aj=exp(fθ(xj))>0a_j= (f_θ(x_j))>0 and Si(q)=∑j:yj≥yiajqjS_i(q)= _j:y_j≥ y_ia_jq_j. The weighted Cox objective is ℒ(q,θ)=∑iδiqilogSi(q)−logaiL(q,θ)= _i _iq_i\ S_i(q)- a_i\, with zero-mass terms interpreted by continuity. Lemma 3 (Sufficient conditions for weighted Cox). Let 0<α<10<α<1, q′∈n−1q _n-1 with qi′>0q _i>0 for all i, and b=q′/αb=q /α. Suppose the weighted Cox objective has a unique global maximizer q⋆q over ℱ(b)F(b), Si(q⋆)>0S_i(q )>0 for every observed event i with δi=1 _i=1, and the coordinates of g⋆=∇qℒ(q⋆,θ)g = _qL(q ,θ) are pairwise distinct. Then the objective is continuous on the simplex and continuously differentiable near q⋆q , and the global inner optimizer satisfies the local stability condition (14). Consequently, the supergradient characterization in Proposition 1 applies. Proof. We first establish continuity of ℒL, continuous differentiability near q⋆q , and convergence of global maximizers as the upper bounds approach b. We then use the strict ordering of the gradient coordinates to derive an explicit representation of these maximizers and establish local Lipschitz stability. Since Si(q)≥aiqiS_i(q)≥ a_iq_i, the term qilogSi(q)q_i S_i(q) tends to zero whenever Si(q)→0S_i(q)→ 0 on the simplex. This gives a continuous extension of the objective to the simplex. Positivity of all event risk-set sums at q⋆q ensures continuous differentiability in a neighborhood of q⋆q . The point q′q is strictly feasible relative to the simplex because 0<qi′<bi0<q _i<b_i. Let d=mini(bi−qi′)>0d= _i(b_i-q _i)>0. For sufficiently small Δb b, the point (1−t)q⋆+tq′(1-t)q +tq , where t=‖Δb‖∞/dt=\| b\|_∞/d, belongs to ℱ(b+Δb)F(b+ b) and converges to q⋆q . By compactness and continuity, any limit of global maximizers for nearby bounds is therefore a global maximizer at b. Uniqueness implies that all such maximizers converge to q⋆q . Choose a permutation π such that gπ(1)⋆>⋯>gπ(n)⋆g _π(1)>·s>g _π(n). Continuity of the gradient ensures that this ordering persists at every nearby global maximizer q q. The first-order necessary condition on the convex set ℱ(b+Δb)F(b+ b) gives ⟨∇qℒ(q^,θ),q−q^⟩≤0,q∈ℱ(b+Δb). _qL( q,θ),q- q ≤ 0, q (b+ b). Thus q q maximizes its linearized objective over this set. With strictly ordered coefficients, this linear problem has the unique solution q^π(k)=minbπ(k)+Δbπ(k),[1−∑j<k(bπ(j)+Δbπ(j))]+,k∈[n]. q_π(k)= \b_π(k)+ b_π(k),\; [1- _j<k(b_π(j)+ b_π(j)) ]_+ \, k∈[n]. This expression is a continuous piecewise affine function of b+Δb+ b, and hence is locally Lipschitz. It equals q⋆q at Δb=0 b=0, proving (14). ∎ We next express the multipliers in terms of the loss gradient. Theorem 2 (Outer supergradient characterization). Fix θ∈Θθ∈ and q′∈n−1q _n-1. Let Φα(q′,θ) _α(q ;θ) be the inner value function defined in (13), where ℒ(⋅,θ):n−1→ℝL(·,θ):S_n-1 is continuous. Assume that the inner maximization problem maxqℒ(q,θ)s.t.∑i=1nqi=1,qi≥0,qi≤1αqi′ _qL(q,θ) .t. _i=1^nq_i=1,\;q_i≥ 0,\;q_i≤ 1αq _i admits an optimal solution q⋆q , and that a constraint qualification holds so that the KKT conditions are satisfied with multipliers ν⋆∈ℝ,λ⋆∈ℝ+n,μ⋆∈ℝ+n,ν , λ ^n_+, μ ^n_+, where ν⋆ν , λ⋆λ , and μ⋆μ are the Lagrange multipliers associated with the constraints ∑i=1nqi=1 _i=1^nq_i=1, qi≥0q_i≥ 0, and qi≤1αqi′q_i≤ 1αq _i, respectively. Assume also that ℒ(⋅,θ)L(·,θ) is continuously differentiable on a neighborhood of q⋆q . Let g⋆:=∇qℒ(q⋆,θ)g := _qL(q ,θ). Then the following statements hold: (i) Stationarity and complementary slackness g⋆−ν⋆+λ⋆−μ⋆=0,g -ν 1+λ -μ =0, λi⋆qi⋆=0,μi⋆(qi⋆−1αqi′)=0,∀i.λ _iq _i=0, μ _i (q _i- 1αq _i )=0, ∀ i. (i) Interior characterization For any i such that 0<qi⋆<1αqi′0<q_i < 1αq _i, ν⋆=gi⋆.ν =g_i . (i) Multiplier characterization For each nondegenerate coordinate with qi′>0q _i>0, μi⋆=(gi⋆−ν⋆)+,if qi⋆=1αqi′,0,if qi⋆<1αqi′.μ _i= cases(g _i-ν )_+,&if q _i= 1αq _i,\\ 0,&if q _i< 1αq _i. cases (iv) Supergradient If, in addition, 0<α<10<α<1, qi′>0q _i>0 for all i, and the local stability condition (14) holds at q⋆=q⋆(b)q =q (b) with b=q′/αb=q /α, then 1αμ⋆∈∂+Φα(q′,θ). 1αμ ∈∂^+ _α(q ;θ). Proof. We first use stationarity and complementary slackness to identify the upper-bound multipliers. Applying Proposition 1 under the additional assumptions in part (iv) then gives the supergradient characterization. Set b:=1αq′,ℱ(b):=q∈ℝn:∑i=1nqi=1,q≥0,q≤b.b:= 1αq , (b):= \q ^n: _i=1^nq_i=1,\ q≥ 0,\ q≤ b \. By the definition of Φα _α, Φα(q′,θ)=maxq∈ℱ(b)ℒ(q,θ). _α(q ;θ)= _q (b)L(q,θ). By assumption, an optimal solution exists and a constraint qualification holds, and hence the KKT conditions are satisfied. Step 1: KKT conditions. Let q⋆q be an optimal solution of the inner problem, and let ν⋆∈ℝ,λ⋆∈ℝ+n,μ⋆∈ℝ+nν , λ ^n_+, μ ^n_+ be the corresponding Lagrange multipliers associated with the constraints ∑i=1nqi=1 _i=1^nq_i=1, qi≥0q_i≥ 0, and qi≤biq_i≤ b_i, respectively. Define g⋆:=∇qℒ(q⋆,θ).g := _qL(q ,θ). Then the stationarity condition is g⋆−ν⋆+λ⋆−μ⋆=0,g -ν 1+λ -μ =0, where ∈ℝn1 ^n denotes the vector of all ones. The complementary slackness conditions are λi⋆qi⋆=0,μi⋆(qi⋆−bi)=0,∀i∈[n]. _i q_i =0, _i (q_i -b_i)=0, ∀ i∈[n]. Step 2: Identification of the multipliers for the upper-bound constraints. Consider an index i such that 0<qi⋆<bi0<q_i <b_i. Then both the lower-bound and upper-bound constraints are inactive, so λi⋆=0,μi⋆=0. _i =0, _i =0. Hence, by the stationarity condition, gi⋆=ν⋆.g_i =ν . If qi⋆<biq_i <b_i, then complementary slackness gives μi⋆=0. _i =0. If 0<qi⋆=bi0<q_i =b_i, then the lower-bound constraint is inactive, so λi⋆=0. _i =0. The stationarity condition then gives μi⋆=gi⋆−ν⋆. _i =g_i -ν . Since μi⋆≥0 _i ≥ 0, we obtain μi⋆=(gi⋆−ν⋆)+,if 0<qi⋆=bi,0,if qi⋆<bi. _i = cases(g_i -ν )_+,&if 0<q_i =b_i,\\ 0,&if q _i<b_i. cases The degenerate case qi⋆=bi=0q_i =b_i=0 is not included in this characterization. In this case, the lower-bound and upper-bound constraints are both active, and the multipliers λi⋆ _i and μi⋆ _i are not separately identifiable from the stationarity condition alone. Step 3: Supergradient characterization. Under the additional assumptions in part (iv), Proposition 1 gives 1αμ⋆∈∂+Φα(q′,θ), 1αμ ∈∂^+ _α(q ;θ), which completes the proof. ∎ Supergradient estimation. At an exact inner maximizer, Theorem 2(i) gives gi⋆=ν⋆g_i =ν on interior coordinates. For an approximate maximizer q q, let g^=∇qℒ(q^,θ) g= _qL( q,θ) and estimate this common value by ν:=1|ℐ|∑i∈ℐg^i,ℐ=i∈[n]:0<q^i<1αqi′.ν:= 1|I| _i g_i, = \i∈[n]:0< q_i< 1αq _i \. Substituting q q, g g, and ν into the multiplier formula gives the supergradient estimate used in the q′q -update. The exact characterization requires the assumptions of part (iv) at an inner optimum. If ℐI is empty, q′q is left unchanged. We use small positive lower bounds on sample weights and weighted risk-set sums for numerical stability. Theorem 3 below establishes local convergence of the outer variables (q′,θ)(q ,θ) under exact inner maximization and the stated assumptions. For computational efficiency, we approximate the inner maximizer using a finite number of gradient steps and estimate the KKT multipliers from the resulting iterate for the outer update. A.1.3 Computational complexity Let n be the sample size, d the covariate dimension, T the number of inner weight updates, and K the number of outer iterations in Algorithm 3. For the linear Cox model, evaluating the risk scores and their parameter derivatives requires O(nd)O(nd) time. Weighted risk-set sums are computed cumulatively after sorting observations by time. The two capped-simplex projections can be computed exactly in O(nlogn)O(n n) time by sorting; our implementation finds the multiplier λ by bisection, which costs O(n)O(n) per bisection step. Consequently, the total computational complexity is O(KTn(d+logn)).O\! (KTn(d+ n) ). For fixed K, T, and d, this cost is nearly linear in the sample size. The additional memory required for the sample weights and cumulative risk-set quantities is O(n)O(n), beyond storage for the data and model parameters. Empirical model-fitting times are reported in Appendix A.3.11. A.1.4 Local convergence analysis We analyze local convergence of the outer variables q′q and θ, assuming that before each q′q -update, we obtain the exact global inner maximizer q⋆q given the current q′q and θ. Under the assumptions below, we show that the outer iterates of the specified alternating scheme with KL-projected q′q -updates and gradient updates of θ converge R-linearly in the Euclidean norm to a strict local minimizer of the resulting outer value function, for a sufficiently small common step size and sufficiently close feasible initialization. Write x=(q′,θ)x=(q ,θ), =1−ϵ(p¯)C=Q_1-ε( p), and F(x)=Φα(q′,θ)F(x)= _α(q ;θ). Let q⋆(x)q (x) denote the global inner maximizer and μ⋆(x)μ (x) its upper-bound KKT multipliers. With DKL(u∥v)=∑iuilog(ui/vi)D_KL(u\|v)= _iu_i (u_i/v_i), the outer weight update takes the entropic mirror-descent form of Beck and Teboulle (2003): q′+=argminu∈⟨μ⋆(x)α,u−q′⟩+1ηDKL(u∥q′).q += _u \ μ (x)α,u-q + 1ηD_KL(u\|q ) \. This is equivalent to the exponentiated update (7) with KL projection. We retain the update order of Algorithm 1: T weight steps with θ fixed, followed by a parameter step using the same inner solution as the last weight step. For the analysis, set ηq′=ηθ=η _q = _θ=η and fix T≥1T≥ 1. Theorem 3 (Local convergence). Let 0<α<10<α<1, 0<ϵ<10<ε<1, and x⋆=(q′⋆,θ⋆)x =(q ,θ ) be feasible, with q′⋆>0q >0 and θ⋆θ interior to the parameter space. Assume: (i) Smooth inner solution. ℒL is twice continuously differentiable near (q⋆(x⋆),θ⋆)(q (x ),θ ). For every x near x⋆x , the global inner maximizer q⋆(x)q (x) is unique, at least one coordinate satisfies 0<qi⋆(x)<qi′/α0<q_i (x)<q _i/α, and both q⋆(x)q (x) and μ⋆(x)μ (x) are continuously differentiable in x. (i) Outer stationarity. The outer KKT conditions hold at x⋆x , and every active upper bound of q′q has a strictly positive outer KKT multiplier. (i) Positive local curvature. At x⋆x , the Hessian of F is positive definite along directions that preserve the simplex constraint and keep the active bound coordinates of q′q fixed, with θ free to vary. Then (q′⋆,θ⋆)(q ,θ ) is a strict local minimizer of the outer value function F. For a sufficiently small fixed η>0η>0 and feasible initialization sufficiently close to (q′⋆,θ⋆)(q ,θ ), there exist constants C>0C>0 and ρ∈(0,1)ρ∈(0,1) such that the outer iterates of the specified KL-projected alternating scheme satisfy ‖(q′(k)−q′⋆θ(k)−θ⋆)‖2≤Cρk,k=0,1,2,…. \| pmatrixq (k)-q \\ θ^(k)-θ pmatrix \|_2≤ Cρ^k, k=0,1,2,…. Proof. We first show that F is smooth and x⋆x is a strict local minimizer. Next, we show that the KL update fixes the weight coordinates whose upper bounds are active at x⋆x . With these coordinates fixed, we prove that a full update cycle contracts the error in a weighted norm, yielding the stated Euclidean bound. Step 1: Smoothness and local minimality. By assumption (i), the inner optimizer varies smoothly with q′q and θ. Differentiating the optimized loss and applying the inner KKT conditions gives, for perturbations along the simplex, dF=⟨μ⋆(x)α,dq′⟩+⟨∇θℒ(q⋆(x),θ),dθ⟩.dF= μ (x)α,dq + _θL(q (x),θ),dθ . The derivatives in this expression are continuously differentiable by assumption (i), so F is twice continuously differentiable locally. Outer stationarity and strict complementarity in (i), together with the positive curvature in (i), give the second-order sufficient conditions for x⋆x to be a strict local minimizer. Step 2: Identification of active upper bounds. We next show that nearby KL updates place the same coordinates of q′q at their upper bounds as at x⋆x . We adapt the active-constraint argument of Bertsekas (1976) to the KL update. At x⋆x , its subproblem has a positive definite objective Hessian, independent active constraint gradients, and strictly positive active multipliers by (i). Its solution therefore varies smoothly with nearby inputs while retaining the same active upper bounds. These weight coordinates can thus be held fixed when analyzing subsequent nearby updates. Step 3: Contraction and the Euclidean error bound. With the active coordinates fixed, we analyze the remaining weight coordinates jointly with θ. Let U have orthonormal columns spanning the directions that preserve ∑iqi′=1 _iq _i=1 and keep qi′q _i fixed whenever qi′⋆=p¯i/(1−ϵ)q_i = p_i/(1-ε). Write q′=q′⋆+Uzq =q +Uz, so z describes changes in the remaining weight coordinates. By assumption (i), the Hessian H of F(q′⋆+Uz,θ)F(q +Uz,θ) at (0,θ⋆)(0,θ ) is positive definite. Set Mq=U⊤diag(1/q′⋆)U,M=diag(Mq,I),D=diag(Mq−1,0).M_q=U diag(1/q )U, M=diag(M_q,I), D=diag(M_q^-1,0). Zero-dimensional weight blocks are omitted. Differentiating the KL optimality equations gives R=I−ηDHR=I-η DH for a weight-only step and J=I−ηM−1HJ=I-η M^-1H for the final weight and parameter step. The derivative of a full cycle is therefore JRT−1JR^T-1. For sufficiently small η>0η>0, H≻0H 0 and 0⪯D⪯M−10 D M^-1 imply ‖J‖H<1\|J\|_H<1 and ‖R‖H≤1\|R\|_H≤ 1, where ‖v‖H2=v⊤Hv\|v\|_H^2=v Hv and the corresponding matrix norm is induced by this vector norm. By continuity of the cycle map’s derivative, the full cycle is therefore a local contraction in this norm after these upper-bound coordinates have been fixed. Iterating the contraction gives geometric decay of the error in the H-norm. Norm equivalence and the orthonormality of U yield the stated Euclidean bound CρkCρ^k, with all fixed factors and the finite initial segment absorbed into C. Continuity of q⋆(x)q (x) also gives convergence of the exact inner maximizers. ∎ A.1.5 Alternating gradient algorithm Algorithm 3 alternates updates of the adversarial weights q and nominal weights q′q , followed by a model-parameter update. Algorithm 3 Alternating gradient algorithm for the coupled DRO problem 1: Input: step sizes ηq,ηq′,ηθ _q, _q , _θ; update steps T, K; initial values θ(0),q(0),q′(0)θ^(0),q^(0),q (0) 2: for k=0,1,2,…,K−1k=0,1,2,…,K-1 do 3: Set q(k,0)←q(k)q^(k,0)← q^(k), q′(k,0)←q′(k)q (k,0)← q (k) 4: ⊳ Update (q,q′)(q,q ) 5: for t=0,…,T−1t=0,…,T-1 do 6: Compute g(k,t)←∇qℒ(q(k,t),θ(k))g^(k,t)← _qL(q^(k,t),θ^(k)) 7: Set q~i(k,t+1)←qi(k,t)exp(ηqgi(k,t))∑j=1nqj(k,t)exp(ηqgj(k,t)) q^(k,t+1)_i← q^(k,t)_i ( _qg^(k,t)_i ) _j=1^nq^(k,t)_j ( _qg^(k,t)_j ) 8: Project q(k,t+1)←Πα(q′(k,t))(q~(k,t+1))q^(k,t+1)← _Q_α(q (k,t))( q^(k,t+1)) 9: Compute g(k,t+1)←∇qℒ(q(k,t+1),θ(k))g^(k,t+1)← _qL(q^(k,t+1),θ^(k)) 10: Define ℐ(k,t+1)←i∈[n]:0<qi(k,t+1)<1αqi′(k,t)I^(k,t+1)← \i∈[n]:0<q^(k,t+1)_i< 1αq (k,t)_i \ 11: if ℐ(k,t+1)=∅I^(k,t+1)= then 12: Set q′(k,t+1)←q′(k,t)q (k,t+1)← q (k,t) 13: else 14: Compute ν(k,t+1)←1|ℐ(k,t+1)|∑i∈ℐ(k,t+1)gi(k,t+1)ν^(k,t+1)← 1|I^(k,t+1)| _i ^(k,t+1)g^(k,t+1)_i 15: Set si(k,t+1)←1α(gi(k,t+1)−ν(k,t+1))+,if 0<qi(k,t+1)=1αqi′(k,t),0,otherwise.s^(k,t+1)_i← cases 1α(g^(k,t+1)_i-ν^(k,t+1))_+,&if 0<q^(k,t+1)_i= 1αq (k,t)_i,\\ 0,&otherwise. cases 16: Set q~i′(k,t+1)←qi′(k,t)exp(−ηq′si(k,t+1))∑j=1nqj′(k,t)exp(−ηq′sj(k,t+1)) q (k,t+1)_i← q (k,t)_i (- _q s^(k,t+1)_i ) _j=1^nq (k,t)_j (- _q s^(k,t+1)_j ) 17: Project q′(k,t+1)←Π1−ϵ(p¯)(q~′(k,t+1))q (k,t+1)← _Q_1-ε( p)( q (k,t+1)) 18: end if 19: end for 20: Set q(k+1)←q(k,T)q^(k+1)← q^(k,T), q′(k+1)←q′(k,T)q (k+1)← q (k,T) 21: ⊳ Update θ 22: Update θ(k+1)←θ(k)−ηθ∇θℒ(q(k+1),θ(k))θ^(k+1)←θ^(k)- _θ _θL(q^(k+1),θ^(k)) 23: end for A.2 Sensitivity analysis A.2.1 Robustness parameters We examine sensitivity to α and ϵε, which control robustness to subpopulation shift and outlier contamination, respectively. We consider α∈0.05,0.2,0.4,1α∈\0.05,0.2,0.4,1\ and ϵ∈0,0.01,0.05,0.1ε∈\0,0.01,0.05,0.1\. When varying α, we fix ϵ=0ε=0; when varying ϵε, we fix α=0.1α=0.1. We use the data-generating process in Section 6.1 with adjusted mixture proportions and Ntrain=Ntest=500N_train=N_test=500. The training distribution is Ptrain=0.95P1+0.05P2P_ train=0.95P_1+0.05P_2 for subpopulation shift and P~train=0.7P1+0.2P2+0.1P~ P_ train=0.7P_1+0.2P_2+0.1 P for outlier contamination. Both settings use the clean test distribution Ptest=0.8P1+0.2P2P_ test=0.8P_1+0.2P_2. Figure 2: Sensitivity to α and ϵε. (Top-left and top-right) Overall and worst-group test C-index for varying α under subpopulation shift. (Bottom) Overall test C-index for varying ϵε under outlier contamination. In Figure 2, smaller α improves worst-group C-index, with a slight decrease in overall C-index at α=0.05α=0.05. Under contamination, positive ϵε reduces training oscillations, with smoother trajectories at ϵ=0.05ε=0.05 and 0.10.1. A.3 Additional experiments and experimental details We examine robustness in additional simulated and real-data settings and describe the experimental procedures. All experiments were conducted on a server equipped with Intel Xeon Gold 6230R CPUs and NVIDIA GeForce RTX 3090 GPUs. A.3.1 Simulation: Subpopulation shift We isolate subpopulation shift by using uncontaminated mixtures of two latent groups. The training distribution is Ptrain=π1trP1+π2trP2P_ train=π^tr_1P_1+π^tr_2P_2, with (π1tr,π2tr)=(0.95,0.05)(π^tr_1,π^tr_2)=(0.95,0.05) and Ntrain=3,000N_train=3,000. The test distribution is Ptest=π1teP1+π2teP2P_ test=π^te_1P_1+π^te_2P_2, with (π1te,π2te)=(0.8,0.2)(π^te_1,π^te_2)=(0.8,0.2) and Ntest=5,000N_test=5,000. The remaining data-generation settings follow Section 6.1. Figure 3: Overall and worst-group C-index under subpopulation shift for the Cox model, f-DRO (left panels), and the proposed method (right panels), across robustness parameters. Box plots summarize 20 independent runs. Figure 3 shows that both robust methods can improve worst-group C-index at a cost in overall C-index, depending on the robustness parameters. A.3.2 Simulation: Nonlinear risk functions We extend the simulation in Section 6.1 to nonlinear group-specific risk functions with linear, quadratic, and interaction terms. We index the majority, minority, and outlier components by g=0,1,2g=0,1,2, respectively, and define r0(x) r_0(x) =x1+x2+0.15x12+0.15x22+0.15x1x2, =x_1+x_2+0.15x_1^2+0.15x_2^2+0.15x_1x_2, r1(x) r_1(x) =x1−x2+0.15x12−0.15x22−0.15x1x2, =x_1-x_2+0.15x_1^2-0.15x_2^2-0.15x_1x_2, r2(x) r_2(x) =−3x1−3x2−0.45x12−0.45x22−0.45x1x2. =-3x_1-3x_2-0.45x_1^2-0.45x_2^2-0.45x_1x_2. Survival times follow the group-specific hazard hg(t∣x)=h0exp(rg(x))h_g(t x)=h_0 (r_g(x)). The training mixture proportions are (0.90,0.08,0.02)(0.90,0.08,0.02), and the clean test proportions are (0.80,0.20,0)(0.80,0.20,0). We use a single-hidden-layer neural Cox model with width 5; the training and test sample sizes, covariate generation, and censoring setup are otherwise the same as in Section 6.1. Training trajectories. Figure 4 summarizes the test C-index trajectories, including the substantial variation of f-DRO and the proposed method with ϵ=0ε=0. With ϵ=0.05ε=0.05, the proposed method attains higher late-stage worst-group C-index than Cox-N and f-DRO, with lower overall C-index than Cox-N. The corresponding IBS trajectories are provided in Appendix A.3.4. Figure 4: Test C-index under nonlinear risks and outlier contamination. Points show interval means, whiskers show minima and maxima over epochs, and lines connect means. Epoch 70 is shown separately. Repeated-run performance. Table 4 reports a separate experiment averaged over five independent repetitions. The proposed method with ϵ=0.05ε=0.05 outperforms f-DRO on all four metrics and improves worst-group C-index and IBS relative to Cox-N, at a cost in overall performance. Table 4: Test performance under nonlinear risks and outlier contamination, averaged over five independent repetitions. Bold and underlined entries denote the best and second-best results, respectively. Method Overall C-index ↑ Worst-group C-index ↑ Overall IBS ↓ Worst-group IBS ↓ Cox-N 0.753 0.562 0.125 0.241 f-DRO 0.551 0.522 0.161 0.217 Ours (ϵ=0.05ε=0.05) 0.698 0.686 0.152 0.199 A.3.3 Simulation: Higher dimension, smaller samples, and heavier censoring We evaluate the methods with d=10d=10 covariates, training sample sizes n∈500,1000n∈\500,1000\, 5,000 test samples, and 50% censoring. Without contamination, the training mixture proportions of the two clean groups are (0.9,0.1)(0.9,0.1); with contamination, the majority, minority, and outlier proportions are (0.90,0.08,0.02)(0.90,0.08,0.02). The clean test proportions remain (0.80,0.20)(0.80,0.20) in both settings. We define a shared risk component and a component controlling group heterogeneity as C(x)=0.6x1+0.5x5,S(x)=1.2x2+0.72(x3+x4)+0.35x6.C(x)=0.6x_1+0.5x_5, S(x)=1.2x_2+ 0.7 2(x_3+x_4)+0.35x_6. The majority and minority log-risk functions are r0(x)=C(x)+0.7S(x),r1(x)=C(x)−0.7S(x),r_0(x)=C(x)+0.7S(x), r_1(x)=C(x)-0.7S(x), and the outlier log-risk function is r2(x)=−2(C(x)+0.7S(x))r_2(x)=-2(C(x)+0.7S(x)). The covariates x7,…,x10x_7,…,x_10 are noise variables and do not affect the true risk. The robust subpopulation parameters are fixed at δ=α=0.1δ=α=0.1 for f-DRO and the proposed method, respectively. For the proposed method, we compare ϵ∈0,0.01,0.02,0.05ε∈\0,0.01,0.02,0.05\ under contamination. Tables 5 and 6 report overall and worst-group C-index, IBS, and tail loss. Without contamination, the robust methods improve worst-group C-index and IBS relative to Cox, with lower overall C-index and higher overall IBS. Under contamination, increasing ϵε improves overall C-index and IBS, while intermediate values achieve the best worst-group C-index and IBS among the evaluated settings. Table 5: Simulation results with d=10d=10, n=500n=500 training samples, and 50% censoring. Values are mean (standard deviation) over 10 runs. Bold and underlined entries denote the best and second-best results within each contamination setting, respectively. Method (n=500n=500) Overall C-index ↑ Worst-group C-index ↑ Overall IBS ↓ Worst-group IBS ↓ Overall tail loss ↓ Worst-group tail loss ↓ Without Outliers Cox 0.787 (0.007) 0.458 (0.025) 0.150 (0.004) 0.357 (0.022) 11.315 (0.373) 11.850 (0.581) f-DRO 0.716 (0.043) 0.537 (0.063) 0.194 (0.016) 0.230 (0.015) 8.711 (0.130) 8.309 (0.059) Ours (ϵ=0ε=0) 0.708 (0.038) 0.526 (0.049) 0.207 (0.012) 0.230 (0.008) 8.568 (0.068) 8.245 (0.019) With Outliers Cox 0.787 (0.005) 0.456 (0.022) 0.150 (0.003) 0.351 (0.012) 10.957 (0.164) 11.298 (0.286) f-DRO 0.628 (0.044) 0.469 (0.088) 0.212 (0.012) 0.246 (0.016) 9.030 (0.148) 8.694 (0.244) Ours (ϵ=0ε=0) 0.599 (0.045) 0.469 (0.104) 0.224 (0.002) 0.242 (0.008) 8.844 (0.154) 8.397 (0.025) Ours (ϵ=0.01ε=0.01) 0.732 (0.019) 0.475 (0.051) 0.189 (0.005) 0.248 (0.009) 8.858 (0.121) 8.340 (0.059) Ours (ϵ=0.02ε=0.02) 0.737 (0.038) 0.522 (0.061) 0.182 (0.014) 0.237 (0.024) 8.819 (0.096) 8.359 (0.078) Ours (ϵ=0.05ε=0.05) 0.764 (0.015) 0.458 (0.050) 0.169 (0.010) 0.287 (0.026) 9.275 (0.248) 8.574 (0.153) Table 6: Simulation results with d=10d=10, n=1000n=1000 training samples, and 50% censoring. Values are mean (standard deviation) over 10 runs. Bold and underlined entries denote the best and second-best results within each contamination setting, respectively. Method (n=1000n=1000) Overall C-index ↑ Worst-group C-index ↑ Overall IBS ↓ Worst-group IBS ↓ Overall tail loss ↓ Worst-group tail loss ↓ Without Outliers Cox 0.787 (0.007) 0.454 (0.017) 0.148 (0.005) 0.344 (0.013) 10.877 (0.169) 11.169 (0.233) f-DRO 0.701 (0.061) 0.523 (0.081) 0.210 (0.010) 0.228 (0.011) 8.722 (0.102) 8.425 (0.135) Ours (ϵ=0ε=0) 0.693 (0.047) 0.525 (0.089) 0.217 (0.005) 0.227 (0.013) 8.535 (0.075) 8.264 (0.025) With Outliers Cox 0.782 (0.008) 0.449 (0.024) 0.153 (0.008) 0.340 (0.028) 10.464 (0.368) 10.468 (0.577) f-DRO 0.557 (0.091) 0.422 (0.077) 0.224 (0.019) 0.246 (0.017) 8.827 (0.168) 8.545 (0.227) Ours (ϵ=0ε=0) 0.592 (0.104) 0.441 (0.083) 0.223 (0.008) 0.238 (0.008) 8.503 (0.011) 8.270 (0.021) Ours (ϵ=0.01ε=0.01) 0.716 (0.050) 0.522 (0.060) 0.202 (0.008) 0.227 (0.008) 8.593 (0.013) 8.249 (0.078) Ours (ϵ=0.02ε=0.02) 0.746 (0.025) 0.455 (0.063) 0.188 (0.008) 0.244 (0.016) 8.755 (0.093) 8.224 (0.017) Ours (ϵ=0.05ε=0.05) 0.774 (0.009) 0.460 (0.045) 0.166 (0.006) 0.271 (0.020) 9.135 (0.126) 8.317 (0.119) A.3.4 Simulation: Integrated Brier score trajectories Figures 5 and 6 report the integrated Brier score (IBS) across training for the linear and nonlinear simulations, respectively. The settings follow Section 6.1 and Appendix A.3.2. Worst-group IBS is the maximum IBS across the two clean test groups. At epoch 70, Cox and Cox-N achieve the lowest overall IBS in their respective settings, while the proposed method with ϵ=0.05ε=0.05 achieves the lowest worst-group IBS. Figure 5: Test IBS under linear risks and outlier contamination. Points show interval means, whiskers show minima and maxima over epochs, and lines connect means. Epoch 70 is shown separately. Figure 6: Test IBS under nonlinear risks and outlier contamination, using the plotting conventions in Figure 5. A.3.5 Benchmark datasets and contamination setup We describe the benchmark populations, data splits, and contamination settings used in the real-data experiments. METABRIC The METABRIC dataset contains breast cancer survival data with genomic, transcriptomic, and clinical information (Curtis et al., 2012; Liu et al., 2014). We use 13 clinical covariates. The dataset contains 1,7931,793 samples, among which 764764 observations are censored, resulting in a censoring rate of 42.61%42.61\%. We define four groups by estrogen receptor (ER) status and age at a threshold of 6060 years. Group 00 contains ER-negative patients with age ≤60≤ 60, Group 11 contains ER-negative patients with age >60>60, Group 22 contains ER-positive patients with age ≤60≤ 60, and Group 33 contains ER-positive patients with age >60>60. The group sizes are 281281, 134134, 525525, and 853853, respectively. The corresponding censoring rates are 48.40%48.40\%, 35.82%35.82\%, 62.67%62.67\%, and 29.43%29.43\%. FLC The FLC dataset (Dispenzieri et al., 2012) comes from a community cohort of Olmsted County, Minnesota residents aged 5050 years or older and was used to study serum free light chains and mortality. We use 66 covariates following Hu and Chen (2024). The dataset contains 7,8747,874 samples, among which 5,7055,705 observations are censored, resulting in a censoring rate of 72.45%72.45\%. We define four groups by sex and age at a threshold of 6060 years. Group 00 consists of females with age ≤60≤ 60, Group 11 consists of males with age ≤60≤ 60, Group 22 consists of females with age >60>60, and Group 33 consists of males with age >60>60. The group sizes are 1,7681,768, 1,6251,625, 2,5822,582, and 1,8991,899, respectively. The corresponding censoring rates are 92.82%92.82\%, 90.09%90.09\%, 59.80%59.80\%, and 55.61%55.61\%. Data split Each dataset is split into training and test sets at an 80:20 ratio, stratified by group labels and censoring indicators. The training set is then split into training and validation subsets at the same ratio, stratified only by censoring indicators. Group labels are used only for the initial split and worst-group evaluation, not for training or model selection. Outlier construction Outliers are introduced only into the training data; validation and test sets remain clean. For linear models, we use a risk-based adversarial construction. We first train a standard Cox model on the training data and compute sample-wise risk scores ri=fθ^Cox(xi).r_i=f_ θ_Cox(x_i). (15) We then select the top 2% highest-risk training samples and construct perturbed risks by sign inversion, r~i=−ri. r_i=-r_i. (16) We generate replacement survival times from the perturbed hazard model T~i∼P~(⋅∣xi),h~(t∣xi)=h0(t)exp(r~i). T_i P(· x_i), h(t x_i)=h_0(t) ( r_i). (17) The generated times replace the selected samples’ original times, creating outcomes that conflict with the fitted risk ordering. The Cox model is fitted with an elastic-net penalty (0.010.01, mixing ratio 0.50.5) on standardized covariates, so the risk scores are centered; the ranking includes censored and uncensored samples, and the event indicators of the selected samples are unchanged. The baseline hazard is constant, h0(t)=λ0h_0(t)= _0 with λ0=∑iδi/∑iyieri _0= _i _i/ _iy_ie^r_i, the maximum-likelihood rate given the fitted scores. For nonlinear models, we use the same risk-based construction with the top 5% highest-risk training samples on both datasets. A.3.6 Baselines We compare methods that require no group labels for training or model selection. • Elastic-Cox: A Cox proportional hazards model with elastic-net regularization. We tune the elastic-net regularization strength over 0,0.1,0.2,0.4,0.8\0,0.1,0.2,0.4,0.8\. • Cox-N: A neural-network Cox model with one hidden layer and hidden dimension 2. We tune the Cox-N regularization strength over 0,0.1,0.2,0.4,0.8\0,0.1,0.2,0.4,0.8\. • f-DRO (Hu and Chen, 2024): A survival DRO method based on the χ2χ^2-divergence and sample splitting. We tune the DRO robustness parameter over 0.1,0.2,0.4,0.8\0.1,0.2,0.4,0.8\. • Exact-Cox (Hu and Chen, 2024): The exact DRO-Cox variant, which reparameterizes the Cox model with a piecewise-constant baseline hazard so that the full likelihood decomposes across observations, and applies the χ2χ^2-DRO objective to the resulting per-observation losses without sample splitting. This construction is specific to the Cox model. We tune the DRO robustness parameter over 0.1,0.2,0.4,0.8\0.1,0.2,0.4,0.8\, initialize the baseline-hazard parameters at their maximum-likelihood values for a zero risk score, and use the same training budget and validation criterion as f-DRO. • Robust-Cox (Minder and Bednarski, 1996): A robust proportional hazards estimator for outlier contamination, applied to linear and nonlinear Cox models. We use the R package coxrobust for linear models. For nonlinear models, we implement an extension of the original estimator by replacing its linear predictor with the same one-hidden-layer network. Starting from a standard neural Cox fit, we optimize the truncated-weight partial likelihood over three refinement rounds. We tune its truncation parameter over 0.90,0.95,0.975,0.99\0.90,0.95,0.975,0.99\ and hold the remaining settings fixed. A.3.7 Training and hyperparameter selection Neural-network training. Cox-N, f-DRO, Exact-Cox, and the proposed method share a one-hidden-layer network of width 2 (ReLU hidden layer, linear output), except in the network configurations evaluated in Table 3, and all are trained with Adam for 100 full-batch parameter updates (K=100K=100 outer iterations for the proposed method) and evaluated at the final iterate. They differ only in the step size (0.010.01 for the baselines, 0.10.1 for the proposed method) and in regularization: Cox-N tunes an elastic-net penalty over 0,0.1,0.2,0.4,0.8\0,0.1,0.2,0.4,0.8\, f-DRO and Exact-Cox use none, and the proposed method fixes it at 0.010.01. For linear models, f-DRO and Exact-Cox run 100 Adam steps with step size 0.050.05, whereas the proposed method runs K=200K=200 gradient steps with step size 0.10.1 and a fixed elastic-net penalty of 0.10.1. The nonlinear simulation uses the architecture in Appendix A.3.2. Optimization details of the proposed method. Algorithm 1 uses T=1T=1 inner weight update per outer iteration and K=200K=200 outer iterations. Step sizes are ηq=0.01 _q=0.01, ηq′=0.1 _q =0.1, and ηθ=0.1 _θ=0.1; the θ-step is plain gradient descent for the linear model and Adam for the nonlinear model. An elastic-net penalty (strength 0.10.1 linear, 0.010.01 nonlinear, mixing ratio 0.50.5) is added to the weighted objective. Projections onto 1−ϵ(p¯)Q_1-ε( p) and α(q′)Q_α(q ) are KL projections onto the capped simplex, [Π(x)]i=minci,λxi[ (x)]_i= \c_i,λ x_i\ with c the coordinate-wise upper bounds and λ found by bisection so that the result sums to one. Validation criteria. Validation and test losses are Cox negative partial log-likelihoods. Elastic-Cox and Cox-N are selected by maximizing validation log-likelihood; Robust-Cox and the DRO methods by minimizing validation tail loss, the average of the largest 5% or 10% of validation losses (Duchi and Namkoong, 2021). The grids for the proposed method are α∈0.1,0.2,0.4,0.6,0.8α∈\0.1,0.2,0.4,0.6,0.8\ and ϵ∈0,0.05,0.1ε∈\0,0.05,0.1\; simulation grids are given in the corresponding subsections. Practical choice of search ranges. Prior knowledge of the smallest meaningful subgroup and of the expected contamination proportion can guide the ranges for α and ϵε; otherwise clustering, anomaly detection, or data-quality diagnostics provide heuristic ranges, and the final pair is selected by validation tail loss. IBS evaluation. We compute IBS from Breslow survival curves with inverse-probability-of-censoring weights estimated by Kaplan–Meier on the training set. Integration uses 80 equally spaced points from the smallest observed time to the 90th percentile of test times, capped at the maximum training time; test samples beyond training follow-up are excluded only from IBS evaluation. A.3.8 Real data analysis: Subpopulation shift Table 7 reports results under subpopulation shift without outliers. Table 7: Results on METABRIC and FLC under subpopulation shift without outliers. The DRO-based methods use the 10% validation tail loss for model selection; Elastic-Cox and Cox-N use validation log-likelihood. Values are mean (standard deviation) over 10 independent runs. METABRIC Model Method IBS ↓ Test log-likelihood ↑ Overall C-index ↑ Worst-group C-index ↑ Test tail loss ↓ Linear Elastic-Cox 0.167 (0.004) -5.003 (0.037) 0.664 (0.014) 0.575 (0.047) 6.432 (0.162) f-DRO 0.183 (0.005) -5.117 (0.033) 0.671 (0.014) 0.583 (0.034) 5.829 (0.017) Exact-Cox 0.173 (0.007) -5.044 (0.064) 0.656 (0.034) 0.550 (0.053) 6.109 (0.085) Ours 0.186 (0.004) -5.128 (0.034) 0.683 (0.015) 0.607 (0.029) 5.829 (0.009) Nonlinear Cox-N 0.167 (0.005) -5.004 (0.043) 0.670 (0.015) 0.584 (0.027) 6.537 (0.103) f-DRO 0.182 (0.005) -5.106 (0.032) 0.666 (0.028) 0.574 (0.041) 5.841 (0.054) Exact-Cox 0.179 (0.006) -5.102 (0.040) 0.627 (0.034) 0.525 (0.083) 6.056 (0.053) Ours 0.179 (0.004) -5.078 (0.028) 0.655 (0.026) 0.575 (0.042) 5.912 (0.034) FLC Model Method IBS ↓ Test log-likelihood ↑ Overall C-index ↑ Worst-group C-index ↑ Test tail loss ↓ Linear Elastic-Cox 0.093 (0.004) -6.467 (0.071) 0.797 (0.004) 0.629 (0.039) 8.087 (0.262) f-DRO 0.111 (0.012) -6.801 (0.200) 0.798 (0.005) 0.629 (0.038) 7.400 (0.091) Exact-Cox 0.115 (0.007) -6.875 (0.112) 0.792 (0.007) 0.620 (0.038) 7.323 (0.065) Ours 0.120 (0.002) -6.956 (0.022) 0.798 (0.004) 0.631 (0.036) 7.283 (0.006) Nonlinear Cox-N 0.092 (0.001) -6.452 (0.028) 0.797 (0.005) 0.608 (0.026) 8.268 (0.036) f-DRO 0.109 (0.011) -6.765 (0.183) 0.797 (0.006) 0.623 (0.032) 7.434 (0.097) Exact-Cox 0.116 (0.003) -6.890 (0.041) 0.774 (0.021) 0.556 (0.078) 7.330 (0.024) Ours 0.114 (0.003) -6.839 (0.048) 0.780 (0.028) 0.605 (0.051) 7.357 (0.061) The proposed method reduces test tail loss relative to the standard Cox models in all settings, while these baselines achieve better IBS and test log-likelihood. It attains the highest overall and worst-group C-index on METABRIC with linear models; on FLC, its worst-group C-index is within 0.0030.003 of the standard Cox models. Rationale for controlled contamination. Luo et al. (2022) model survival outliers through individual changes in hazard rates, allowing survival times to be unusually long or short relative to predictions based on covariates. Our constructions introduce such departures at a known contamination proportion while preserving the covariates and clinically defined evaluation groups; their effects can differ across datasets and model classes. All methods use the same contamination construction within each setting. Appendix A.3.9 evaluates an alternative that appends perturbed observations while preserving the original training samples. A.3.9 Additional outlier setting We evaluate linear models on FLC after appending perturbed copies of short-survival event samples, preserving the original training data. Let ℰ=i:zi∈train,δi=1E=\i:z_i _train,\ _i=1\ (18) denote the index set of uncensored event samples in the training data. We sort the samples indexed by ℰE by their observed survival times yiy_i in ascending order and select the lowest 2%2\%. Specifically, we define low=i∈ℰ:yi≤c0.02.O_low= \i :y_i≤ c_0.02 \. (19) Here, c0.02c_0.02 is the empirical 2nd percentile of the observed event times among samples in ℰE. Each selected sample i∈lowi _low is duplicated with its covariates and event indicator unchanged, and its survival time resampled near the median event time: x~i=xi,δ~i=δi,y~i∼(mE−τ,mE+τ), x_i=x_i, δ_i= _i, y_i (m_E-τ,m_E+τ), (20) where mEm_E is the median observed survival time among event samples and τ=0.1mEτ=0.1m_E. The resulting duplicated sample is denoted by z~i=(x~i,y~i,δ~i). z_i=( x_i, y_i, δ_i). (21) The contaminated training set is ~train=train∪z~i:i∈low. D_train=D_train∪\ z_i:i _low\. (22) The appended samples retain covariates from short-survival cases but have times near the median, creating inconsistencies between covariates and outcomes. Table 8: Linear-model results on FLC with appended outliers. The 10% and 5% settings specify the validation tail loss for the DRO methods; Elastic-Cox uses validation log-likelihood. Values are mean (standard deviation) over 10 independent runs. FLC Model Method IBS ↓ Test log-likelihood ↑ Overall C-index ↑ Worst-group C-index ↑ Test tail loss ↓ Linear (10%) Elastic-Cox 0.091 (0.001) -6.430 (0.025) 0.798 (0.004) 0.618 (0.030) 7.520 (0.017) f-DRO 0.115 (0.016) -7.377 (1.157) 0.728 (0.106) 0.579 (0.092) 8.933 (2.887) Ours 0.113 (0.002) -6.842 (0.031) 0.797 (0.004) 0.644 (0.028) 7.180 (0.009) Linear (5%) Elastic-Cox 0.091 (0.001) -6.430 (0.025) 0.798 (0.004) 0.618 (0.030) 7.520 (0.017) f-DRO 0.129 (0.006) -7.629 (0.011) 0.688 (0.089) 0.563 (0.087) 10.177 (4.925) Ours 0.120 (0.002) -6.970 (0.018) 0.797 (0.005) 0.646 (0.031) 7.283 (0.004) In Table 8, Elastic-Cox achieves the best IBS, test log-likelihood, and overall C-index. The proposed method achieves the best worst-group C-index and test tail loss under both validation criteria, with overall C-index close to Elastic-Cox. For f-DRO, the 5% validation criterion yields lower overall C-index and higher test tail loss than the 10% criterion. A.3.10 Contamination sensitivity Tables 1 and 3 use the same risk-inversion construction (Appendix A.3.5) on both datasets, with 2% contamination for linear models and 5% for nonlinear models. To examine sensitivity to the contamination rate, Table 9 evaluates all four model–dataset combinations at both rates, using the 10% validation setting of Table 1. The (linear, 2%) and (nonlinear, 5%) entries coincide with the 10% validation results in Tables 1 and 3. At 5%, the proposed method attains the highest worst-group C-index in all four settings. At 2%, it leads in the linear settings but trails the best nonlinear baselines. Its test tail loss is the lowest or second lowest in every setting. Table 9: Worst-group C-index and 5% test tail loss under the risk-based contamination applied to all model classes and datasets at contamination rates of 2% and 5%. Hyperparameters are selected by the 10% validation tail loss for Robust-Cox, f-DRO, Exact-Cox, and the proposed method and by validation log-likelihood for Elastic-Cox and Cox-N. Values are mean (standard deviation). 2% contamination 5% contamination Model Method Worst-group C-index ↑ Test tail loss ↓ Worst-group C-index ↑ Test tail loss ↓ METABRIC Linear Elastic-Cox 0.501 (0.058) 6.305 (0.111) 0.466 (0.064) 6.195 (0.178) Robust-Cox 0.562 (0.042) 6.229 (0.048) 0.520 (0.052) 6.117 (0.091) f-DRO 0.574 (0.030) 5.825 (0.015) 0.567 (0.031) 5.822 (0.015) Exact-Cox 0.459 (0.057) 6.179 (0.164) 0.393 (0.032) 6.267 (0.176) Ours 0.613 (0.028) 5.835 (0.011) 0.609 (0.028) 5.835 (0.008) Nonlinear Cox-N 0.572 (0.038) 6.404 (0.084) 0.508 (0.091) 6.561 (0.310) Robust-Cox 0.586 (0.040) 6.351 (0.077) 0.545 (0.049) 6.295 (0.161) f-DRO 0.559 (0.043) 5.832 (0.049) 0.533 (0.047) 5.839 (0.041) Exact-Cox 0.419 (0.105) 6.139 (0.100) 0.395 (0.078) 6.226 (0.178) Ours 0.571 (0.039) 5.901 (0.029) 0.585 (0.040) 5.902 (0.014) FLC Linear Elastic-Cox 0.618 (0.038) 7.650 (0.114) 0.602 (0.035) 7.560 (0.128) Robust-Cox 0.611 (0.040) 7.854 (0.031) 0.594 (0.040) 7.626 (0.030) f-DRO 0.626 (0.037) 7.409 (0.035) 0.618 (0.035) 7.384 (0.048) Exact-Cox 0.609 (0.037) 7.581 (0.127) 0.584 (0.044) 7.506 (0.173) Ours 0.631 (0.036) 7.283 (0.006) 0.632 (0.036) 7.284 (0.006) Nonlinear Cox-N 0.614 (0.037) 7.917 (0.183) 0.623 (0.029) 7.801 (0.208) Robust-Cox 0.626 (0.041) 8.071 (0.173) 0.602 (0.030) 8.134 (0.324) f-DRO 0.626 (0.031) 7.457 (0.035) 0.624 (0.028) 7.447 (0.043) Exact-Cox 0.582 (0.057) 7.603 (0.209) 0.561 (0.081) 7.622 (0.188) Ours 0.613 (0.040) 7.359 (0.064) 0.630 (0.033) 7.371 (0.067) A.3.11 Computational runtime Table 10 compares model-fitting times for linear Cox models on METABRIC and FLC. For this comparison, all optimization-based methods, including the proposed method, are run for 100 iterations, separately from the training schedules in Appendix A.3.7. The proposed method is faster than f-DRO on both datasets, although standard Cox and Exact-Cox are faster. With the K=200K=200 schedule of Appendix A.3.7, it takes 0.6310.631 and 1.2831.283 seconds on METABRIC and FLC, respectively. The corresponding complexity analysis is given in Appendix A.1.3. Table 10: Model-fitting time in seconds, reported as mean (standard deviation) over 10 repetitions on a single CPU thread. All optimization-based methods are run for 100 iterations in this comparison. Dataset Cox Elastic-Cox f-DRO Exact-Cox Ours METABRIC 0.104 (0.005) 1.239 (0.040) 1.129 (0.281) 0.195 (0.004) 0.368 (0.008) FLC 0.356 (0.133) 2.210 (0.796) 18.122 (0.531) 0.236 (0.007) 0.817 (0.135)