Paper deep dive

Probabilistic Multi-Source Data Selection via Learnable Categorical Priors: Variance Reduction and Guaranteed Derandomization

Shizhan Cai, Qixin Zhang, Siyuan Ma, Siyuan Liang, Yongcheng Jing, Dacheng Tao

Year: 2026Venue: arXiv preprintArea: cs.LGType: PreprintEmbeddings: 181

Intelligence

Status: not_run | Model: - | Prompt: - | Confidence: 0%

Entities (0)

No extracted entities yet.

Relation Signals (0)

No relation signals yet.

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Bilevel data selection is a fundamental task with numerous applications in machine learning and statistics. A common strategy for tackling this problem is to relax the original bilevel combinatorial optimization into a continuous one. However, existing methods, including 0-1 relaxation and probabilistic relaxation, still suffer from several limitations, such as unstable gradient estimation and infeasible rounding. To address these issues, this paper proposes a new probabilistic relaxation for the multi-source bilevel data selection problem. Our key idea is to introduce multiple learnable categorical priors to parameterize a probability distribution supported only on feasible data subsets, thus naturally yielding a lossless rounding scheme. We further show that our proposed relaxation admits an effective two-point gradient estimator, which can simultaneously avoid the inverse-Hessian computation required by existing 0-1 relaxation and the unbounded score term appearing in prior probabilistic relaxation. Furthermore, we develop a reducible representation of this estimator, which substantially reduces the number of lower-level retrainings required for gradient computation. Built upon this reducible two-point estimator, we then design a novel stochastic gradient algorithm, termed ReSGD, for multi-source bilevel data selection problem, and establish a rigorous approximation guarantee for the subset returned by ReSGD under weak submodularity. Finally, extensive experiments are conducted to validate the effectiveness of our proposed ReSGD.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

177,292 characters extracted from source content.

Expand or collapse full text

_set:Ne hessianbox _set:Ne _set:Ne Probabilistic Multi-Source Data Selection via Learnable Categorical Priors: Variance Reduction and Guaranteed Derandomization Shizhan Cai Qixin Zhang Siyuan Ma Affiliation: Siyuan Liang Yongcheng Jing Dacheng Tao Affiliation: College of Computing and Data Science Affiliation: Nanyang Technological University Email: shizhan001@e.ntu.edu.sg qixin.zhang2026@gmail.com Email: masi0004@e.ntu.edu.sg Email: siyuan.liang,yongcheng.jing,dacheng.tao@ntu.edu.sg Abstract Bilevel data selection is a fundamental task with numerous applications in machine learning and statistics. A common strategy for tackling this problem is to relax the original bilevel combinatorial optimization into a continuous one. However, existing methods, including 0-1 relaxation and probabilistic relaxation, still suffer from several limitations, such as unstable gradient estimation and infeasible rounding. To address these issues, this paper proposes a new probabilistic relaxation for the multi-source bilevel data selection problem. Our key idea is to introduce multiple learnable categorical priors to parameterize a probability distribution supported only on feasible data subsets, thus naturally yielding a lossless rounding scheme. We further show that our proposed relaxation admits an effective two-point gradient estimator, which can simultaneously avoid the inverse-Hessian computation required by existing 0-1 relaxation and the unbounded score term appearing in prior probabilistic relaxation. Furthermore, we develop a reducible representation of this estimator, which substantially reduces the number of lower-level retrainings required for gradient computation. Built upon this reducible two-point estimator, we then design a novel stochastic gradient algorithm, termed ReSGD, for multi-source bilevel data selection problem, and establish a rigorous approximation guarantee for the subset returned by ReSGD under weak submodularity. Finally, extensive experiments are conducted to validate the effectiveness of our proposed ReSGD. 1 Introduction Data selection aims to identify a small subset of training samples such that training solely on this subset can achieve performance comparable to that trained on the full dataset. As modern learning systems continue to grow in scale, this task has become increasingly crucial in a wide range of applications, including efficient training (Chang et al., 2026; Fan et al., 2025; Fan et al., 2026; Huang et al., 2024; Mindermann et al., 2022; Mirzasoleiman et al., 2020), in-context learning (Kumari et al., 2024a; Kumari et al., 2024b; Zhang et al., 2025) and continual learning (Borsos et al., 2020; Borsos et al., 2024; Hao et al., 2024; Tong et al., 2025; Zhou et al., 2022). For instance, in continual learning, an effective way to counter catastrophic forgetting (French, 1999; McCloskey and Cohen, 1989) is to maintain a small rehearsal memory that stores a subset of previously seen data and then replay them during training on new tasks (Buzzega et al., 2020; Chaudhry et al., 2019; Lopez-Paz and Ranzato, 2017; Yoon et al., 2022). From an optimization viewpoint, many of data selection problems (Borsos et al., 2020; Chen et al., 2025b; Killamsetty et al., 2021; Lin et al., 2024) can be formulated as bilevel combinatorial optimization, where the upper level evaluates the selected subset according to the performance of the model induced by the corresponding lower-level training problem. However, the discrete nature of data selection makes the underlying bilevel problem highly challenging. To circumvent this, recent works have sought to bring the original bilevel combinatorial optimization into the world of continuous optimization, thereby leveraging the rich toolbox of gradient-based methods (Bertsekas, 2015; Lan, 2020; Nesterov, 2013). Existing approaches mainly follow two routes. The first is 0-1 relaxation (Borsos et al., 2020; Hao et al., 2024), which directly relaxes the binary selective mask into a continuous variable in [0,1][0,1] and then applies implicit differentiation (Lorraine et al., 2020; Pedregosa, 2016) to solve the resulting bilevel continuous optimization. The second is probabilistic relaxation (Zhou et al., 2022), which treats each binary mask as a Bernoulli random variable with learnable parameters and then optimizes these parameters via policy gradient descent Williams (1992); Sutton et al. (1999); Kakade (2001). Although these two aforementioned continuous-relaxation reformulations have made data selection tasks more tractable, their practical implementation still faces several fundamental limitations: ➊ Problematic Gradient Estimations. In 0-1 relaxation, estimating the upper-level gradient typically requires approximating the inverse Hessian of the lower-level objective, which is often computationally prohibitive (Borsos et al., 2020; Hao et al., 2024; Zhang et al., 2024). In contrast, probabilistic relaxation avoids such second-order computation via policy-gradient estimators Sutton et al. (1999); Williams (1992), but these estimators are unbounded and often suffer from instability and high variance during training (Greensmith et al., 2004; Wu et al., 2018); ➋ Limited to Single-Source Scenarios. Both 0-1 and probabilistic relaxation primarily focus on subset selection over a single homogeneous dataset. However, in many real-world applications, we often encounter multiple heterogeneous data sources rather than a single-source dataset. A representative example arises in continual learning, where one often needs to retain a small data subset separately from each past task for future replay (Buzzega et al., 2020; Lopez-Paz and Ranzato, 2017; Tong et al., 2025). Such settings cannot be trivially decomposed into multiple independent single-source data selection problems, since this treatment will ignore the cross-source correlations among datasets and lead to redundant selections; ➌ Heuristic and Infeasible Rounding. In continuous relaxation, rounding or derandomization refers to the conversion of the obtained continuous solution into a discrete solution (Ageev and Sviridenko, 2004; Bu et al., 2024). However, satisfactory rounding schemes are still lacking for both 0-1 and probabilistic relaxation. In 0-1 relaxation, current rounding methods are heuristic, so there is no principled way to quantify how the obtained subset relates to the original continuous solution. In probabilistic relaxation, sampling from the learned Bernoulli distributions may return a subset outside the target constraint (Calinescu et al., 2011; Chekuri et al., 2010; Chekuri et al., 2014); ➍ Incomplete Convergence Analysis. Existing theoretical results of both 0-1 and probabilistic relaxation mainly focus on convergence to a continuous stationary point (Zhou et al., 2022; Hao et al., 2024; Chen et al., 2025b; Borsos et al., 2024). However, this is insufficient for data selection, since our primary interest lies in the quality of the resulting subset rather than the intermediate continuous solution. To date, how to establish rigorous theoretical guarantees for the obtained subset remains an uncharted territory in bilevel data selection. In view of all this, a natural question arises: Can we develop a provable continuous-relaxation framework, together with guaranteed rounding scheme and effective gradient estimator, for multi-source data selection problems? In this paper, we will provide an affirmative answer to this question by presenting a novel probabilistic relaxation for multi-source data selection problems. Our key idea is to introduce multiple learnable categorical priors to parameterize a probability distribution supported only on feasible data subsets. Compared with prior relaxations, this constrained construction immediately yields a lossless rounding scheme for multi-source data selection problem. Moreover, we derive an effective two-point gradient estimator for our proposed probabilistic relaxation, which can simultaneously avoid the costly inverse-Hessian computation in 0-1 relaxation (Borsos et al., 2020; Hao et al., 2024) and the unbounded score in prior probabilistic relaxation (Zhou et al., 2022), thereby achieving a significant variance-reduction effect. Based on this two-point estimator, we further develop an SGD algorithm, termed ReSGD, for the multi-source data selection problem, and establish a rigorous approximation guarantee for the subset returned by ReSGD under weak submodularity (Das and Kempe, 2018; Thiery and Ward, 2022). In summary, we make the following contributions: ∙ We present a novel probabilistic relaxation for multi-source data selection problem, which only assigns probability mass to feasible subsets and thus naturally induces a lossless rounding scheme. ∙ We prove that our proposed probabilistic relaxation admits an effect two-point gradient estimator, which avoids both the inverse-Hessian computation in 0-1 relaxation and the unbounded score in prior probabilistic relaxation, thereby improving stability and reducing variance. Furthermore, we develop a reducible representation of this two-point gradient estimator under a Bayesian perspective, which can further reduce the gradient-computation cost from ⁡(h+∑i=1hni)O(h+ _i=1^hn_i) retrainings to ⁡(h)O(h). Here, h denotes the number of data sources, and ∑i=1hni _i=1^hn_i is the total sample size. ∙ Based on the proposed reducible two-point representation, we then design a novel SGD algorithm, termed ReSGD. Furthermore, under (γ,β)(γ,β)-weak submodularity (Thiery and Ward, 2022), we show that the subset returned by ReSGD achieves a (γ2β+β⁡(1−γ)+γ2) ( γ^2β+β(1-γ)+γ^2 )-approximation guarantee with respect to the optimal data subset ∗S^*. To the best of our knowledge, this provides the first rigorous theoretical guarantee for multi-source bilevel data selection problem. ∙ We demonstrate the practical efficacy of our ReSGD on data summarization and continual learning. 2 Bilevel Data Selection: From Single Dataset to Multi-Source Scenarios This section first reviews the single-source bilevel data selection problem together with its two representative continuous relaxations, namely, 0-1 relaxation (Borsos et al., 2020; Borsos et al., 2024; Hao et al., 2024) and probabilistic relaxation (Zhou et al., 2022). We then move on to the multi-source setting, which is the main focus of this paper. 2.1 Single-Source Bilevel Data Selection: Background and Relaxations Given a fixed dataset ≜(xi,yi)i=1nD \(x_i,y_i)\_i=1^n, the goal of single-source data selection is to identify a subset of at most K samples from D such that training on this subset can achieve performance comparable to that trained on D. Let ℓ⁡(xi,yi,) (x_i,y_i; θ) be the training loss on the i-th sample (xi,yi)(x_i,y_i), where θ denotes the model parameters. Then, the single-source data selection problem can be formulated as the following bilevel combinatorial optimization (Borsos et al., 2020; Borsos et al., 2024; Dumouchelle et al., 2024; Hao et al., 2024; Tong et al., 2025; Zhou et al., 2024; Zhou et al., 2022): min∈0,1n,‖0≤K _m∈\0,1\^n,~\|m\|_0≤ K ℒ⁡(∗())≜∑i=1nℓ⁡(xi,yi,∗()) ( θ^*(m)) _i=1^n (x_i,y_i; θ^*(m)) (1) s.t.∗()∈ .t. θ^*(m)∈ arg⁡minℒ^(,)≜∑i=1nmi⋅ℓ⁡(xi,yi,), _ θ~ L( θ,m) _i=1^nm_i· (x_i,y_i; θ), where ∈0,1nm∈\0,1\^n denotes the selective mask vector with mi=1m_i=1 indicating that the i-th sample is selected and mi=0m_i=0 otherwise, ∗() θ^*(m) is the optimal solution to the lower-level problem under the selective mask vector m and the symbol ∥⋅∥0\|·\|_0 represents the ℓ0 _0 norm. The bilevel problem in Eq. (1) is generally known to be NP-hard (Feige, 1998; Natarajan, 1995), implying that no polynomial-time algorithms can solve it optimally. To tackle this, recent works have relaxed Eq. (1) into a continuous one, so that gradient-based methods can be used to approximately solve it. To date, two representative relaxations have been developed for Eq. (1), namely, 0-1 and probabilistic relaxation. 0-1 Relaxation. The first line of work (Borsos et al., 2020; Borsos et al., 2024; Hao et al., 2024) directly relaxes each mask mim_i as a continuous variable over [0,1][0,1]. In this way, Eq. (1) can be rewritten as a bilevel continuous optimization, that is, min∈[0,1]n,‖1≤K _m∈[0,1]^n,~\|m\|_1≤ K ℒ⁡(∗())≜∑i=1nℓ⁡(xi,yi,∗()), ( θ^*(m)) _i=1^n (x_i,y_i; θ^*(m)), (2) s.t.∗()∈ .t. θ^*(m)∈ arg⁡minℒ^(,)≜∑i=1nmi⋅ℓ⁡(xi,yi,). _ θ~ L( θ,m) _i=1^nm_i· (x_i,y_i; θ). Compared with Eq. (1), Eq. (2) enables the use of well-established gradient-based methods from bilevel continuous optimization (Chen et al., 2025a; Hong et al., 2023; Ji et al., 2021) to search for a high-quality fractional mask vector m. In particular, under suitable conditions, the hypergradient w.r.t. each mi∈[0,1]m_i∈[0,1] can be written as: ∂ℒ⁡(∗())∂mi≜−∇ℒ(∗())⊤ [∂2ℒ^(∗(),)∂∂T]−1∇ℓ(xi,yi,∗()). ( θ^*(m))∂ m_i - _ θL( θ^*(m)) [ ∂^2 L( θ^*(m),m)∂ θ∂ θ^T ]^-1 _ θ (x_i,y_i; θ^*(m)). (3) It is worth noting that Eq. (3) requires computing the inverse Hessian matrix of the lower-level objective function ℒ L, which often incurs substantial time and memory overhead. Probabilistic Relaxation. To avoid the costly inverse Hessian computation in Eq. (3), (Zhou et al., 2022) recently proposed a probabilistic relaxation, which models each mask mim_i as an independent Bernoulli random variable, namely, mi∼Bern⁡(πi)m_i ( _i) where πi∈[0,1] _i∈[0,1] is a learnable parameter. Let ≜(π1,…,πn) π ( _1,…, _n) denote the corresponding parameter vector. Under the independence assumption, the induced distribution over binary vector m is p(∣)≜∏i=1nπimi(1−πi)1−mip(m π) _i=1^n _i^m_i(1- _i)^1-m_i. Accordingly, Eq. (1) can be relaxed into an expected loss minimization problem with respect to the continuous parameter vector ≜(π1,…,πn) π ( _1,…, _n), that is to say, min∈[0,1]n,‖1≤K⁡Φ⁡()≜p(|)[ℒ⁡(∗())]≜∑∈0,1n(ℒ⁡(∗())⋅p(|)) _ π∈[0,1]^n,\| π\|_1≤ K ( π) _p( m| π) [L ( θ^* (m ) ) ] _m∈\0,1\^n (L( θ^*(m))·p( m| π) ) (4) s.t.∗()∈arg⁡min⁡ℒ^(,)≜∑i=1nmi⋅ℓ⁡(xi,yi,). .t. θ^*(m)∈ _ θ\ L( θ,m) _i=1^nm_i· (x_i,y_i; θ). A key advantage of (4) is that its gradient admits an unbiased policy gradient estimator Williams (1992); Sutton et al. (1999),i.e., ∂Φ⁡()∂πi≜p(|)[ℒ⁡(∗())⋅∇πiln(p(|))]≜p(|)[ℒ⁡(∗())⋅ mi-πiπi(1-πi)]. ∂ ( π)∂ _i _p(m| π) [L( θ^*(m))· _ _i (p(m| π)) ] _p(m| π) [L( θ^*(m))· $ m_i- _i _i(1- _i)$ ]. (5) Eq.(5) implies that ℒ⁡(∗())⋅mi−πiπi(1−πi)L( θ^*(m))· m_i- _i _i(1- _i) is an unbiased stochastic gradient estimator of the relaxed objective function Φ⁡() ( π). In sharp contrast with Eq.(3), this estimator only requires evaluating the function value ℒ(∗())L( θ^*(m)) and is thus easier to implement. However, it is worth noting that when mi=1,πi→0m_i=1, _i→ 0 or mi=0,πi→1m_i=0, _i→ 1, the factor mi−πiπi(1−πi) m_i- _i _i(1- _i) will become unbounded, which can lead to high variance, unstable updates, and even divergence in practice (Greensmith et al., 2004; Wu et al., 2018). Furthermore, the Bernoulli product distribution p(∣)p(m π) in Eq. (4) poses a significant challenge for rounding, since it assigns non-zero probability mass to subsets outside the K-size cardinality constraint. As a consequence, for an arbitrary parameter vector π, sampling according to p(∣)p(m π) may yield an infeasible subset whose cardinality exceeds K. 2.2 Multi-Source Bilevel Data Selection We now turn to the multi-source scenarios. Specifically, we first suppose that the full dataset D consists of h heterogeneous sources, that is to say, ≜∪i∈[h]iD _i∈[h]D_i, where each source is given by i≜(xj(i),yj(i))∣j∈[ni]D_i \(x^(i)_j,y^(i)_j) j∈[n_i]\ and i∩j=∅,∀i≠jD_i _j= ,∀ i≠ j. Here, different sources may correspond to different tasks, domains or data providers, depending on the applications. Given a set of predefined budgets K1,…,Kh\K_1,…,K_h\, the goal of multi-source data selection is to select at most KiK_i from each source iD_i such that the model trained on the jointly selected subset can achieve comparable performance as the model trained on D. Formally, we consider the following problem: mini∈0,1ni,‖i‖0≤Ki,∀i∈[h]⁡ℒ⁡(∗())≜∑i=1h∑j=1niℓ⁡(xj(i),yj(i),∗()) _m_i∈\0,1\^n_i,\|m_i\|_0≤ K_i,∀ i∈[h]L( θ^*(m)) _i=1^h _j=1^n_i (x^(i)_j,y^(i)_j; θ^*(m)) (6) s.t.∗()∈arg⁡min⁡ℒ^(,)≜∑i=1h∑j=1nimi,j⋅ℓ⁡(xj(i),yj(i),), .t. θ^*(m)∈ _ θ\ L( θ,m) _i=1^h _j=1^n_im_i,j· (x^(i)_j,y^(i)_j; θ), where the variable mi,j∈0,1m_i,j∈\0,1\ indicates whether the j-th sample from the i-th source iD_i is selected, i≜(mi,1,…,mi,ni)m_i (m_i,1,…,m_i,n_i) and ≜(1,…,h)m (m_1,…,m_h). Previously, when facing the multi-source data selection problem, many works (Borsos et al., 2020; Borsos et al., 2024; Hao et al., 2024; Zhou et al., 2022) artificially decomposed it into multiple independent single-source selection problems and then repeatedly applied cardinality-constrained algorithms to each source. In contrast, the upper-level objective in our Eq. (6) evaluates the model trained on the jointly selected subset across all sources, so the benefit of selecting a sample from one source depends crucially on what has been selected in other sources. 3 A New Probabilistic Relaxation via Learnable Categorical Priors This section mainly introduces a new probabilistic relaxation for the multi-source data selection problem (6). As discussed in Section 2.1, a critical drawback of the current probabilistic relaxation (Zhou et al., 2022) is that the induced Bernoulli product distribution will assign non-zero probability to subsets outside the target constraint, which renders the subsequent rounding operation intrinsically infeasible. To overcome this, a natural idea is to design a new sampling distribution whose support is entirely restricted to the feasible region. In other words, for any combinatorial constraint C, we seek to construct a parameterized distribution p⁡(∣)p(m π) such that p⁡(∣)≡0p(m π)≡ 0 for any ∉m . An immediate advantage of such a construction is that every sample drawn from p⁡(∣)p(m π) is automatically feasible. Next, we show how to instantiate this idea for the multi-source data selection problem (6). Recall that Eq. (6) requires selecting at most KiK_i samples from each source iD_i. From a high-level perspective, this process closely resembles an election procedure, where each source iD_i is viewed as a local candidate pool and KiK_i samples with the highest vote counts are chosen. Motivated by this observation, we introduce a learnable categorical prior i≜(πi,1,…,πi,ni) π_i ( _i,1,…, _i,n_i) for every source i≜(xj(i),yj(i))∣j∈[ni]D_i \(x^(i)_j,y^(i)_j) j∈[n_i]\. Here, each πi,j _i,j represents the probability of selecting the j-th sample from source iD_i, while the residual probability 1−∑j=1niπi,j1- _j=1^n_i _i,j corresponds to selecting a dummy symbol ∅ , meaning that no sample from iD_i is chosen. With these categorical priors ≜(1,…,h) π ( π_1,…, π_h), we then independently draw KiK_i samples from each source iD_i according to the categorical distribution specified by i π_i. This operation will yield a random subset ⋃i=1h⋃k=1Kidi,k _i=1^h _k=1^K_id_i,k, where di,kd_i,k denotes the k-th random draw from source iD_i. Note that this subset ⋃i=1h⋃k=1Kidi,k _i=1^h _k=1^K_id_i,k contains at most KiK_i samples from each source iD_i. Therefore, this sampling rule induces a probability distribution supported only in the feasible set of Eq. (6). Based on this constrained sampling scheme, we naturally can relax Eq. (6) into the following continuous optimization problem over the categorical priors π, that is to say, mini∈Δni,∀i∈[h]⁡Φ⁡() _ π_i∈ _n_i,\,∀ i∈[h] ( π) ≜∑di,k∈i∪∅,∀i,k(ℒ(∗(∪i=1h∪k=1Kidi,k))⋅∏i=1h∏k=1Kip(di,k∣i)), _d_i,k _i∪\ \,∀ i,k (L ( θ^* ( _i=1^h _k=1^K_id_i,k ) )· _i=1^h _k=1^K_ip(d_i,k π_i) ), (7) s.t.∗() .t. θ^*(S) ∈argminℒ^(,)≜∑(x,y)∈ℓ(x,y;),∀⊆, ∈ _ θ~ L( θ,S) _(x,y) (x,y; θ),\,\,\, , where di,kd_i,k denotes the k-th random draw from source iD_i with p⁡((xj(i),yj(i))∣i)≜πi,jp((x^(i)_j,y^(i)_j) π_i) _i,j, Δni≜∈[0,1]ni:‖1≤1 _n_i \ p∈[0,1]^n_i:\ \| p\|_1≤ 1\ denotes the sub-probability simplex and p⁡(∅∣i)≜1−∑j=1niπi,jp( π_i) 1- _j=1^n_i _i,j. Remark 1. Note that, in Eq. (7), we switch from the binary notation used in Eq. (6) to an equivalent set-based formulation. Specifically, for any subset ⊆S , let ∗() θ^*(S) denote the optimal lower-level solution trained on the data subset S. Equivalently, ∗() θ^*(S) corresponds to ∗() θ^*(1_S) of Eq. (6). Compared with the continuous relaxation in Eq. (4), a major advantage of Eq. (7) is that it only assigns probability mass to feasible subsets. Consequently, for any given categorical priors ≜(1,…,h) π ( π_1,…, π_h), we can directly sample a subset satisfying the multi-source budget constraints. Moreover, this sampling step is lossless in expectation w.r.t. the relaxed objective value Φ⁡() ( π). Specifically, if the subset S is sampled via categorical priors π, then [ℒ⁡(∗())]=Φ⁡()E_S [L( θ^*(S)) ]= ( π). 3.1 Variance-Reduced Two-point Gradient Estimator and Reducible SGD To leverage gradient-based methods for optimizing Eq. (7), this subsection focuses on deriving stochastic estimators for ∇Φ() _ π ( π). A natural starting point is to follow the same route as in Eq. (5) and compute a policy gradient estimator for our new probabilistic relaxation Φ⁡() ( π). However, as shown in Appendix C.1, such an estimator still involves an unbounded score factor analogous to that in Eq. (5), which can again lead to high variance and unstable training behavior, especially when the categorical priors approach the boundary of the feasible domain. Fortunately, the special structure of our proposed new probabilistic relaxation Φ⁡() ( π) admits a variance-reduced two-point gradient estimator, which can avoid the aforementioned unbounded score factor. Before going into the details, we first introduce some notations. For any realization of the random draws di,k,∀k∈[Ki],∀i∈[h]\d_i,k,\,∀ k∈[K_i],\,∀ i∈[h]\, we define −(i,1)≜⋃(i′,k′)≠(i,1)di′,k′S_-(i,1) _(i ,k )≠(i,1)d_i ,k as the subset obtained by excluding the first draw di,1d_i,1 from source iD_i. Equivalently, −(i,1)S_-(i,1) is formed by taking the union of all sampled draws except di,1d_i,1. With this notation, we can establish the following theorem: Theorem 1 (Proof in Appendix C.2). For any i∈[h]i∈[h] and j∈[ni]j∈[n_i], the partial derivative of Φ⁡() ( π) with respect to variable πi,j _i,j satisfies the following equality, namely, ∂Φ⁡()∂πi,j≜Ki⋅⁡[ ℒ(∗(−(i,1)∪(xj(i),yj(i))))−ℒ⁡(∗(−(i,1)))⏟Two-point Representation], ∂ ( π)∂ _i,j K_i·E [ ( θ^* (S_-(i,1)∪(x^(i)_j,y^(i)_j) ) )-L ( θ^* (S_-(i,1) ) )_ Two-point Representation ], (8) where the expectation is taken over all random draws except the first draw di,1d_i,1 from source iD_i. Theorem 1 suggests the partial derivative of Φ⁡() ( π) w.r.t. πi,j _i,j can be estimated via the two-point difference ℒ⁡(∗(−(i,1)∪(xj(i),yj(i))))−ℒ⁡(∗(−(i,1)))L( θ^*(S_-(i,1)∪$x^(i)_j,y^(i)_j)$)-L( θ^*(S_-(i,1))). Compared with the previous Eq.(3) and Eq.(5), this two-point representation enjoys three major advantages. First, it only depends on function evaluations of the upper-level objective, thereby avoiding costly backpropagation and the computation of the inverse Hessian matrix required in Eq. (3). Second, compared to Eq.(5), Eq. (8) does not rely on any unbounded factor, thereby leading to more stable updates. Third, the two subsets −(i,1)∪(xj(i),yj(i))S_-(i,1)∪(x^(i)_j,y^(i)_j) and −(i,1)S_-(i,1) only differ in one data point, so their corresponding function values are typically very close to each other. This strong coupling can significantly reduce the fluctuation of the resulting gradient estimator, and thus naturally yields a variance-reduction effect (Duchi et al., 2015; Nesterov and Spokoiny, 2017; Shamir, 2017). It is worth noting that the two-point representation in Eq. (8) requires evaluating the marginal contribution of every data point (xj(i),yj(i))(x_j^(i),y_j^(i)) with respect to its corresponding subset −(i,1)S_-(i,1), which will eventually lead to a total of ⁡(h+∑i=1hni)O(h+ _i=1^hn_i) lower-level retrainings. Here, ∑i=1hni _i=1^hn_i denotes the total number of data points across all sources and is typically very large in practice. As a result, directly computing stochastic gradients based on Eq. (8) is generally impractical. To address this issue, we next introduce a computationally-efficient reformulation of Eq. (8). Algorithm 1 SGD via Two-point Reducible Representation (ReSGD) Input: Datasets ii=1h\D_i\_i=1^h, Budgets Kii=1h\K_i\_i=1^h, Learning rate η, Number of iterations T Output: Selected subset ≜⋃i=1h⋃k=1Kidi,kS _i=1^h _k=1^K_id_i,k Initialize i(0)∈Δni π_i^(0)∈ _n_i for all i∈[h]i∈[h] and train a holdout model ∗() θ^*(D); for t=0,1,…,T−1t=0,1,…,T-1 do Draw KiK_i samples from source iD_i according to i(t) π_i^(t) for all i∈[h]i∈[h]; Construct −(i,1)(t)S^(t)_-(i,1) by excluding the first draw from source iD_i for all i∈[h]i∈[h]; Train ∗(−(i,1)) θ^*(S_-(i,1)) for all i∈[h]i∈[h]; Compute ^i(t)≜(g^i,1(t),…,g^i,ni(t)) g^(t)_i ( g^(t)_i,1,…, g^(t)_i,n_i) by reducible representation Eq. (10) ∀i∈[h],j∈[ni]∀ i∈[h],j∈[n_i]; Update i(t+1)≜Π11footnotemark: 1 Δni(i(t)−η⋅^i(t)),∀i∈[h] π_i^(t+1) to0.0pt 0.80002pt -1.07639pt$^$ _ _n_i ( π_i^(t)-η· g_i^(t) ),∀ i∈[h]; Sample an index τ uniformly at random from 0,1,…,T−1\0,1,…,T-1\; Draw†footnotemark: KiK_i samples di,1,…,di,Ki\d_i,1,…,d_i,K_i\ from source iD_i via i(τ),∀i∈[h] π_i^(τ),∀ i∈[h]; †footnotetext: We provide the details of projection operation and a more effective reservoir rounding scheme in Appendix E and F. Reducible Representation. We now derive a tractable expression for the two-point difference of Eq. (8) under Bayesian statistical theory (Gelman et al., 1995; Murphy, 2012). Specifically, we can treat the model parameter θ as a random variable with prior p⁡()p( θ), and seek to infer its corresponding posterior p⁡(|)p( θ|S) from the already-seen data subset ⊆S . We further suppose that the model follows a predictive distribution p⁡(y|x;)≜∫p⁡(y|x;)⋅p⁡(|)p(y|x;S) (y|x; θ)·p( θ|S)d θ. Note that when a point estimate of θ is used, the posterior p⁡(|)p( θ|S) can be viewed as a Dirac measure (Durrett, 2019; Evans, 2025). Under this probabilistic framework, the loss function is commonly specified as the negative log-likelihood of the predictive distribution, namely, ℓ⁡(x,y,∗())≜−log⁡p⁡(y|x;) (x,y; θ^*(S)) - (y|x;S). Then, throughout the Bayes’ rule and conditional independence, numerous prior works (Mindermann et al., 2022; Tong et al., 2025) have shown that the marginal contribution of any data point (x,y)(x,y) with respect to a subset S can admit the following reducible form:(See Appendix C.3.) ℒ⁡(∗(∪(x,y)))−ℒ⁡(∗())=log⁡p⁡(|;)−log⁡p⁡(|;∪(x,y)) ( θ^* (S∪(x,y) ) )-L ( θ^* (S ) )= ( y| x;S)- ( y| x;S∪(x,y)) (9) =log⁡p⁡(|;)−log⁡p⁡(y|x;)⋅p⁡(|;)p⁡(y|x;)=ℓ⁡(x,y,∗())−ℓ⁡(x,y,∗()), = ( y| x;S)- p(y|x;D)·p( y| x;S)p(y|x;S)= (x,y; θ^*(D) )- (x,y; θ^*(S) ), where the pair (,)( x, y) denotes all samples in D and the second equality comes from Eq.(2) in (Mindermann et al., 2022). By substituting Eq. (9) into Eq. (8), we then can obtain the following reducible representation: ∂Φ⁡()∂πi,j≜Ki⋅⁡[ ℓ(xj(i),yj(i),∗())−ℓ⁡(xj(i),yj(i),∗(−(i,1)))⏟Reducible Two-point Representation]. ∂ ( π)∂ _i,j K_i·E [ (x^(i)_j,y^(i)_j; θ^*(D) )- (x^(i)_j,y^(i)_j; θ^*(S_-(i,1)) )_Reducible Two-point Representation ]. (10) It is worth noting that the reducible two-point representation in Eq. (10) only requires ⁡(h)O(h) lower-level trainings to compute the optimal parameters ⁡(−(i,1)),∀i∈[h] θ(S_-(i,1)),∀ i∈[h] and ⁡() θ(D). This is significantly smaller than the ⁡(h+∑i=1hni)O(h+ _i=1^hn_i) retrainings required by Eq. (8). Furthermore, we emphasize that the above Bayesian probabilistic framework is well motivated, since several prior works have rigorously shown that many neural networks can be modeled via Bayesian probabilistic model (Wild et al., 2023; Liu et al., 2019; Jospin et al., 2022; Gal and Ghahramani, 2016). 4 Stochastic Gradient Descent via Reducible Two-point Representation Based on the reducible two-point representation in Eq. (10), we now present a novel stochastic gradient descent algorithm, termed ReSGD, for multi-source bilevel data selection problem. As summarized in Algorithm 1, ReSGD first initializes a categorical prior i(0)∈Δni π_i^(0)∈ _n_i for each source iD_i and trains a holdout model on the full dataset D. After that, at every iteration t, ReSGD draws KiK_i samples from each source iD_i according to the current categorical prior i(t) π_i^(t). Then, for any i∈[h]i∈[h], we construct −(i,1)(t)S_-(i,1)^(t) by excluding the first draw of iD_i, and train the corresponding model ∗(−(i,1)(t)) θ^*(S_-(i,1)^(t)). With these trained models, ReSGD computes the stochastic gradient vector ^i(t)≜(g^i,1(t),…,g^i,ni(t)) g_i^(t) ( g_i,1^(t),…, g_i,n_i^(t)) via Eq. (10), where each coordinate is given by g^i,j(t)≜Ki[ℓ⁡(xj(i),yj(i),∗())−ℓ⁡(xj(i),yj(i),∗(−(i,1)(t)))] g_i,j^(t) K_i [ (x_j^(i),y_j^(i); θ^*(D) )- (x_j^(i),y_j^(i); θ^*(S_-(i,1)^(t)) ) ]. After obtaining ^i(t) g_i^(t), we update each categorical prior by projected gradient descent. Note that the symbol ΠΔni(⋅) _ _n_i(·) in Algorithm 1 denotes the Euclidean projection onto Δni _n_i. Finally, after running the above SGD iterations for T steps, ReSGD draws KiK_i data samples from each source iD_i and returns their union. 4.1 Approximation Guarantee under Weak Submodularity This subsection aims to theoretically verify the effectiveness of our proposed ReSGD. Specifically, we will exploit an intrinsic combinatorial structure of our concerned multi-source data selection problem, known as weak submodularity (Chen et al., 2018; Das and Kempe, 2018; Elenberg et al., 2018; Thiery and Ward, 2022), to establish an approximation guarantee for ReSGD. Before delving into the details, we first introduce some commonly used notations. Specifically, for any set function f:2→ℝf:2^D , we denote the marginal gain of adding ℬB to A by f⁡(ℬ|)≜f⁡(∪ℬ)−f⁡()f(B|A) f(A )-f(A). Especially when ℬ≜eB \e\ is a singleton, we simply write f⁡(e∣)≜f⁡(e∣)f(e ) f(\e\ ). In addition, we say that a set function f is monotone if and only if f⁡()≤f⁡(ℬ)f(A)≤ f(B) for any ⊆ℬ⊆A . Since weak submodularity is typically used to characterize utility maximization problems, in the remainder of this subsection, we consider an equivalent utility function U(∗())U( θ^*(S)) of the upper-level loss ℒ(∗())L( θ^*(S)) in Eq. (6). Specifically, we define U⁡(∗())≜ℒ⁡(init)−ℒ⁡(∗())U( θ^*(S)) ( θ_init)-L( θ^*(S)), which measures the reduction in the upper-level loss brought by training on S w.r.t. the initial parameter setup init θ_init. With the above notations, we now formally introduce the concept of weak submodularity. Definition 1 ((γ,β)(γ,β)-Weak Submodularity). A utility set function f:2→ℝf:2^D is called (γ,β)(γ,β)-weakly submodular if there exist γ∈(0,1]γ∈(0,1] and β≥1β≥ 1 such that, for any ⊆ℬ⊆A , ∑e∈ℬ∖f⁡(e∣) _e f(e ) ≥γ⋅f(ℬ∣), ≥γ· f(B ), ∑e∈ℬ∖f⁡(e∣ℬ∖e) _e f(e \e\) ≤β⋅f⁡(ℬ∣). ≤β· f(B ). (11) Here, γ is called the lower submodularity ratio, while β is called the upper submodularity ratio. Intuitively, from the perspective of data selection, Eq. (11) measures how closely the accumulated marginal gains of individual data points approximate the joint utility gain of selecting the entire subset. Furthermore, weak submodularity has been widely observed in a broad range of applications, including variable selection (Das and Kempe, 2018; Elenberg et al., 2018), video summarization (Chen et al., 2018; Nguyen and Thai, 2022), neural network pruning (El Halabi et al., 2022; Tofigh et al., 2025), and sparse optimal transport (Manupriya et al., 2024; Alvarez-Melis et al., 2018). Next, we will show that such a combinatorial property also broadly appears in bilevel data selection problems. Warm-up: Bayesian Optimal Design (Chaloner and Verdinelli, 1995; Hashemi et al., 2019; Rainforth et al., 2024). We consider the data selection problem for Bayesian linear regression. Specifically, we assume that the data is generated from a linear model ≜X+ϵ y X θ+ ε, where ≜[y1,…,yn]⊤∈ℝn y [y_1,…,y_n] ^n, X≜[x1,…,xn]⊤∈ℝn×pX [x_1,…,x_n] ^n× p, the noise vector ϵ ε follows a Gaussian distribution ⁡(,σ2In)N( 0,σ^2I_n) and the model parameter θ is equipped with a Gaussian prior, i.e., ∼⁡(,Λ) θ ( 0, ). In this case, the lower-level optimization of bilevel data selection problem can be written as ∗()≜arg⁡min⁡∑(xi,yi)∈(yi−xi⊤)2+σ2⊤Λ−1 θ^*(S) _ θ _(x_i,y_i) (y_i-x_i θ)^2+σ^2 θ ^-1 θ, and our goal is to select a small subset ⊆S such that the resulting posterior estimator has a small prediction error over the full dataset. Equivalently, we seek to minimize this expected upper-level loss function ℒ⁡(∗())≜12nϵ,[‖−X∗()‖22]L( θ^*(S)) 12nE_ ε, θ [\| y-X θ^*(S)\|_2^2 ]. Note that (Borsos et al., 2020; Borsos et al., 2024) recently showed that ℒ(∗())L( θ^*(S)) is closely related to the Bayesian V-optimal criterion ℒV(∗())≜12nϵ,[‖X⁡(−∗())‖22]L_V( θ^*(S)) 12nE_ ε, θ[\|X ( θ- θ^*(S) )\|_2^2]. Specifically, when the number of data points n is sufficiently large, these two objectives differ only by an additive constant σ2/2σ^2/2. Therefore, minimizing ℒ(∗())L( θ^*(S)) can be approximately reduced to minimizing ℒV(∗())L_V( θ^*(S)). More importantly, we further observe that the utility function induced by the Bayesian V-optimal criterion ℒV(∗())L_V( θ^*(S)) is (γ,β)(γ,β)-weakly submodular, as stated below. Theorem 2 (Proof in Appendix B). The utility function UV()≜ℒV()−ℒV(∗())U_V(S) _V( 0)-L_V( θ^*(S)) induced by the Bayesian V-optimal criterion ℒV(∗())L_V( θ^*(S)) is monotone and (1/c,c)(1/c,c)-weakly submodular, where c≜1+s2σ2⋅λmax(Λ)c 1+ s^2σ^2· _max( ), s≜maxi∈[n]⁡‖xi‖2s _i∈[n]\|x_i\|_2 and λmax(Λ) _ ( ) is the largest eigenvalue of the matrix Λ . Motivated by the results of Theorem 2, we next establish an approximation guarantee for our proposed ReSGD algorithm under the weak submodularity. Theorem 3 (Proof in Appendix D). Let n≜||n |D|, Kmax≜maxi∈[h]⁡KiK_max _i∈[h]K_i, learning rate η≜1/nKmax2⋅Tη 1/nK_max^2· T and ∗S^* denotes the optimal solution of Eq. (6). If the utility function U⁡(∗())≜ℒ⁡(init)−ℒ⁡(∗())U( θ^*(S)) ( θ_init)-L( θ^*(S)) induced by Problem Eq. (6) is (γ,β)(γ,β)-weakly monotone submodular, then the subset S returned by Algorithm 1 satisfies ⁡[U⁡(∗())]≥(γ2β+β⁡(1−γ)+γ2)⋅U⁡(∗(∗))−⁡(n⋅Kmax2T).E\! [U( θ^*(S)) ]≥ ( γ^2β+β(1-γ)+γ^2 )· U( θ^*(S^*))-O ( n· K_max^2T ). Remark 2. Theorem 3 implies that, when the induced utility function U(∗())U( θ^*(S)) is (γ,β)(γ,β)-weakly monotone submodular, if we set T≜⁡(nKmax2/ϵ2)T (nK_ ^2/ε^2), the subset yielded by our proposed Algorithm 1 can secure a utility value of (ρ⋅U⁡(∗(∗))−ϵ)(ρ· U( θ^*(S^*))-ε) where ρ≜γ2β+β⁡(1−γ)+γ2ρ γ^2β+β(1-γ)+γ^2. Equivalently, in terms of the original loss function, we have ⁡[ℒ⁡(∗())]≤ρ⋅ℒ⁡(∗(∗))+(1−ρ)⋅ℒ⁡(init)+ϵE\! [L( θ^*(S)) ]≤ρ·L( θ^*(S^*))+(1-ρ)·L( θ_init)+ε. 5 Experiments This section empirically validates the effectiveness of our proposed ReSGD algorithm via two different applications: Data summarization (Borsos et al., 2024; Mirzasoleiman et al., 2013; Schreiber et al., 2020) and Replay-based continual learning (Borsos et al., 2020; Hao et al., 2024; Tong et al., 2025; Zhou et al., 2022). Figure 1: Performance comparison between ReSGD and other baselines on data summarization, illustrating the superior performance and stability of ReSGD under various data corruption scenarios. Data Summarization. We conduct experiments on three benchmarks: MNIST Deng (2012), CIFAR-10 Krizhevsky et al. (2009), and CIFAR-100 Krizhevsky et al. (2009). To ensure a fair evaluation, we follow the model architecture and settings of Borsos et al. (2020); Tong et al. (2025); Zhou et al. (2022). We compare our ReSGD against a set of competitive baselines, including Greedy Coreset Borsos et al. (2020), PBCS Zhou et al. (2022), and CSRel Tong et al. (2025). Additionally, uniform sampling is included to provide a lower-bound performance reference. The experimental details are provided in the Appendix K. Furthermore, to showcase the effectiveness of ReSGD in more practical and challenging scenarios, we follow the data corruptions via label noise (pairwise and symmetric) and class imbalance in Zhou et al. (2022). The details are deferred to Appendix J. As shown in Figure 1, ReSGD consistently outperforms all baselines in terms of test accuracy across different coreset sizes, while exhibiting stronger robustness to data corruptions, including label noise and class imbalance. Similar trends are observed on CIFAR-100, with detailed results provided in Appendix G. Continual Learning. We evaluate ReSGD on Split MNIST Zenke et al. (2017), Perm MNIST Goodfellow et al. (2013), Split CIFAR-10, and Split CIFAR-100. To ensure a fair evaluation, we follow the model architecture and settings of Borsos et al. (2020); Tong et al. (2025); Zhou et al. (2022). We set the memory size to 100 for Split MNIST and Perm MNIST, and 200 for Split CIFAR-10 and Split CIFAR-100. For the Split CIFAR-100 dataset, we split a total of 100 classes into 10 disjoint tasks and use ResNet-18 He et al. (2016) as the backbone. Detailed experiment settings and hyperparameters are shown in Appendix K. The baselines include Uniform, k-Means Features Nguyen et al. (2017), k-Center Embedding Sener and Savarese (2017), Hardest Aljundi et al. (2019), iCaRL Rebuffi et al. (2017), OCS Yoon et al. (2022), BCSR Hao et al. (2024), Greedy Coreset Borsos et al. (2020), PBCS Zhou et al. (2022), and CSReL-CL-Prv Tong et al. (2025). The average accuracies in Table 1 show that our method performs on par with Greedy Coreset on the MNIST datasets and consistently outperforms all other baselines on the more complex datasets (Split CIFAR-10 and Split CIFAR-100), highlighting the superiority of our approach. Additionally, the experimental results indicate that considering previous tasks can further enhance CL performance. Further demonstrations of the effectiveness and compatibility of our method with other CL approaches are deferred to Appendix H. Further Analysis. We conduct effectiveness study in Appendix I and meanwhile provide a detailed sensitivity analysis in Appendix K. Table 1: Performance comparison on continual learning benchmarks. We mark the Best results in bold. Our memory construction method performs on par with Greedy Coreset on the MNIST dataset and consistently outperforms other baselines on more complex datasets. Method PermMNIST SplitMNIST SplitCIFAR-10 SplitCIFAR-100 Uniform 91.82 ± 0.94 78.09 ± 0.80 33.66 ± 1.11 7.06 ± 0.40 k-Means Features 93.22 ± 1.47 78.29 ± 0.59 34.38 ± 1.26 8.05 ± 0.54 k-Center Embedding 93.69 ± 0.59 78.70 ± 0.56 34.07 ± 2.14 8.97 ± 0.47 Hardest 88.17 ± 1.44 77.04 ± 0.84 25.35 ± 1.48 6.43 ± 0.37 iCaRL 93.68 ± 0.45 79.62 ± 0.50 33.26 ± 0.82 8.31 ± 0.22 OCS 84.86 ± 2.69 75.07 ± 0.91 34.12 ± 1.94 9.27 ± 0.45 BCSR 93.81 ± 0.91 78.30 ± 0.81 35.14 ± 3.64 12.11 ± 1.24 Greedy Coreset 95.67 ± 0.21 79.31 ± 0.57 34.88 ± 2.76 8.41 ± 0.25 PBCS 94.22 ± 0.61 76.30 ± 0.81 34.07 ± 2.14 13.02 ± 0.29 CSReL-CL-Prv 95.02 ± 0.23 78.64 ± 0.36 36.82 ± 0.95 15.49 ± 0.84 ReSGD 95.72 ± 0.12 80.03 ± 0.17 37.68 ± 0.13 16.70 ± 0.12 6 Conclusion In this paper, we present a novel probabilistic relaxation for the multi-source data selection problem, which only assigns probability mass to feasible subsets and thus naturally induces a lossless rounding scheme. Most importantly, we show that this relaxation admits an effective two-point gradient estimator, avoiding both the inverse-Hessian computation in 0-1 relaxation and the unbounded score in prior probabilistic relaxation. From a Bayesian perspective, we further derive a reducible representation that reduces the gradient-computation cost from ⁡(h+∑i=1hni)O(h+ _i=1^hn_i) retrainings to ⁡(h)O(h), where h is the number of sources and ∑i=1hni _i=1^hn_i is the total sample size. Building upon this reducible representation, we develop a stochastic algorithm, termed ReSGD and prove that, under (γ,β)(γ,β)-weak submodularity, its returned subset achieves a (γ2β+β⁡(1−γ)+γ2) ( γ^2β+β(1-γ)+γ^2 )-approximation guarantee with respect to the optimal subset ∗S^*. Finally, extensive experiments on data summarization and continual learning validate the effectiveness of our ReSGD. References [1] A. A. Ageev and M. I. Sviridenko (2004) Pipage rounding: a new method of constructing algorithms with proven performance guarantee. Journal of Combinatorial Optimization 8, p. 307–328. Cited by: item ➌. [2] R. Aljundi, K. Kelchtermans, and T. Tuytelaars (2019) Task-free continual learning. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 11254–11263. Cited by: §5. [3] D. Alvarez-Melis, T. Jaakkola, and S. Jegelka (2018) Structured optimal transport. In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, Vol. 84, p. 1771–1780. Cited by: §4.1. [4] A. Attia, S. Leyffer, and T. S. Munson (2022) Stochastic learning approach for binary optimization: application to bayesian optimal design of experiments. SIAM Journal on Scientific Computing 44 (2), p. B395–B427. Cited by: §C.1. [5] A. Attia, M. Schliserman, U. Sherman, and T. Koren (2026) Fast last-iterate convergence of SGD in the smooth interpolation regime. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, External Links: Link Cited by: Appendix L. [6] A. S. Benjamin, D. Rolnick, and K. Kording (2018) Measuring and regularizing networks in function space. arXiv preprint arXiv:1805.08289. Cited by: Appendix H. [7] D. Bertsekas (2015) Convex optimization algorithms. Athena Scientific. Cited by: §1. [8] Z. Borsos, M. Mutny, and A. Krause (2020) Coresets via bilevel optimization for continual learning and streaming. Advances in neural information processing systems 33, p. 14879–14890. Cited by: Appendix A, Figure 5, Figure 5, 1st item, 1st item, §K.1, §K.1, §K.2, §K.2, item ➊, §1, §1, §1, §2.1, §2.1, §2.2, §2, §4.1, §5, §5, §5. [9] Z. Borsos, M. Mutnỳ, M. Tagliasacchi, and A. Krause (2024) Data summarization via bilevel optimization. Journal of Machine Learning Research 25 (73), p. 1–53. Cited by: Appendix A, item ➍, §1, §2.1, §2.1, §2.2, §2, §4.1, §5. [10] J. Bracken and J. T. McGill (1973) Mathematical programs with optimization problems in the constraints. Operations research 21 (1), p. 37–44. Cited by: Appendix A. [11] F. Bu, H. Jo, S. Y. Lee, S. Ahn, and K. Shin (2024) Tackling prevalent conditions in unsupervised combinatorial optimization: cardinality, minimum, covering, and more. In Forty-first International Conference on Machine Learning, External Links: Link Cited by: item ➌. [12] P. Buzzega, M. Boschini, A. Porrello, D. Abati, and S. Calderara (2020) Dark experience for general continual learning: a strong, simple baseline. Advances in neural information processing systems 33, p. 15920–15930. Cited by: Appendix A, Appendix H, item ➋, §1. [13] G. Calinescu, C. Chekuri, M. Pal, and J. Vondrák (2011) Maximizing a monotone submodular function subject to a matroid constraint. SIAM Journal on Computing 40 (6), p. 1740–1766. Cited by: item ➌. [14] K. Chaloner and I. Verdinelli (1995) Bayesian experimental design: a review. Statistical science, p. 273–304. Cited by: §B.1, §4.1. [15] P. Chang, J. Zhang, B. Chen, C. Wang, C. Guo, Y. Zhang, Y. Gao, J. Xiang, Y. Gao, C. Sun, Y. Chen, and D. kong (2026) SPICE: submodular penalized information–conflict selection for efficient large language model training. In The Fourteenth International Conference on Learning Representations, External Links: Link Cited by: §1. [16] A. Chaudhry, P. K. Dokania, T. Ajanthan, and P. H. Torr (2018) Riemannian walk for incremental learning: understanding forgetting and intransigence. In Proceedings of the European conference on computer vision (ECCV), p. 532–547. Cited by: Appendix A. [17] A. Chaudhry, M. Ranzato, M. Rohrbach, and M. Elhoseiny (2018) Efficient lifelong learning with a-gem. arXiv preprint arXiv:1812.00420. Cited by: Appendix H. [18] A. Chaudhry, M. Rohrbach, M. Elhoseiny, T. Ajanthan, P. Dokania, P. Torr, and M. Ranzato (2019) Continual learning with tiny episodic memories. In Workshop on Multi-Task and Lifelong Reinforcement Learning, Cited by: §1. [19] C. Chekuri, J. Vondrak, and R. Zenklusen (2010) Dependent randomized rounding via exchange properties of combinatorial structures. In 2010 IEEE 51st Annual Symposium on Foundations of Computer Science, Vol. , p. 575–584. External Links: Document Cited by: item ➌. [20] C. Chekuri, J. Vondrák, and R. Zenklusen (2014) Submodular function maximization via the multilinear relaxation and contention resolution schemes. SIAM Journal on Computing 43 (6), p. 1831–1879. Cited by: item ➌. [21] L. Chen, Y. Ma, and J. Zhang (2025) Near-optimal nonconvex-strongly-convex bilevel optimization with fully first-order oracles. Journal of Machine Learning Research 26 (109), p. 1–56. Cited by: Appendix A, §2.1. [22] L. Chen, M. Feldman, and A. Karbasi (2018) Weakly submodular maximization beyond cardinality constraints: does randomization help greedy?. In International Conference on Machine Learning, p. 804–813. Cited by: §4.1, §4.1. [23] Z. Chen, Y. Zheng, Z. Lai, Z. Yang, C. Li, Y. Liu, and L. Lin (2025) Quadratic coreset selection: certifying and reconciling sequence and token mining for efficient instruction tuning. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, Cited by: item ➍, §1. [24] A. Das and D. Kempe (2018) Approximate submodularity and its applications: subset selection, sparse approximation and dictionary selection. Journal of Machine Learning Research 19 (3), p. 1–34. Cited by: Appendix A, §1, §4.1, §4.1. [25] L. Deng (2012) The mnist database of handwritten digit images for machine learning research. IEEE Signal Processing Magazine 29 (6), p. 141–142. Cited by: §5. [26] J. C. Duchi, M. I. Jordan, M. J. Wainwright, and A. Wibisono (2015) Optimal rates for zero-order convex optimization: the power of two function evaluations. IEEE Transactions on Information Theory 61 (5), p. 2788–2806. Cited by: §3.1. [27] J. Dumouchelle, E. Julien, J. Kurtz, and E. B. Khalil (2024) Neur2BiLO: neural bilevel optimization. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, External Links: Link Cited by: §2.1. [28] R. Durrett (2019) Probability: theory and examples. Vol. 49, Cambridge university press. Cited by: §3.1. [29] M. El Halabi, S. Srinivas, and S. Lacoste-Julien (2022) Data-efficient structured pruning via submodular optimization. Advances in Neural Information Processing Systems 35, p. 36613–36626. Cited by: §4.1. [30] E. R. Elenberg, R. Khanna, A. G. Dimakis, and S. Negahban (2018) Restricted strong convexity implies weak submodularity. The Annals of Statistics 46 (6B), p. 3539–3568. Cited by: §4.1, §4.1. [31] T. Elsken, J. H. Metzen, and F. Hutter (2019) Neural architecture search: a survey. Journal of Machine Learning Research 20 (55), p. 1–21. Cited by: Appendix A. [32] L. C. Evans (2025) Measure theory and fine properties of functions. Chapman and Hall/CRC. Cited by: §3.1. [33] Z. Fan, S. Du, S. Hu, P. Wang, L. Shen, Z. Ya, D. Tao, and Y. Wang (2025) Combatting dimensional collapse in llm pre-training data via diversified file selection. In The Thirteenth International Conference on Learning Representations, Cited by: §1. [34] Z. Fan, Y. Xian, Y. Sun, K. Shen, and L. Shen (2026) Joint selection for large-scale pre-training data via policy gradient-based mask learning. In The Fourteenth International Conference on Learning Representations, External Links: Link Cited by: §1. [35] U. Feige (1998) A threshold of ln n for approximating set cover. Journal of the ACM 45 (4), p. 634–652. Cited by: §2.1. [36] C. Finn, P. Abbeel, and S. Levine (2017) Model-agnostic meta-learning for fast adaptation of deep networks. In International conference on machine learning, p. 1126–1135. Cited by: Appendix A. [37] L. Franceschi, M. Donini, P. Frasconi, and M. Pontil (2017) Forward and reverse gradient-based hyperparameter optimization. In International conference on machine learning, p. 1165–1173. Cited by: Appendix A. [38] L. Franceschi, P. Frasconi, S. Salzo, R. Grazzi, and M. Pontil (2018) Bilevel programming for hyperparameter optimization and meta-learning. In International conference on machine learning, p. 1568–1577. Cited by: Appendix A. [39] R. M. French (1999) Catastrophic forgetting in connectionist networks. Trends in cognitive sciences 3 (4), p. 128–135. Cited by: Appendix A, §1. [40] Y. Gal and Z. Ghahramani (2016) Dropout as a bayesian approximation: representing model uncertainty in deep learning. In international conference on machine learning, p. 1050–1059. Cited by: §3.1. [41] A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin (1995) Bayesian data analysis. Chapman and Hall/CRC. Cited by: §3.1. [42] S. Ghadimi and M. Wang (2018) Approximation methods for bilevel programming. arXiv preprint arXiv:1802.02246. Cited by: Appendix A. [43] I. J. Goodfellow, M. Mirza, D. Xiao, A. Courville, and Y. Bengio (2013) An empirical investigation of catastrophic forgetting in gradient-based neural networks. arXiv preprint arXiv:1312.6211. Cited by: §K.2, §5. [44] S. Gould, B. Fernando, A. Cherian, P. Anderson, R. S. Cruz, and E. Guo (2016) On differentiating parameterized argmin and argmax problems with application to bi-level optimization. arXiv preprint arXiv:1607.05447. Cited by: Appendix A. [45] R. Grazzi, L. Franceschi, M. Pontil, and S. Salzo (2020) On the iteration complexity of hypergradient computation. In International Conference on Machine Learning, p. 3748–3758. Cited by: Appendix A. [46] E. Greensmith, P. L. Bartlett, and J. Baxter (2004) Variance reduction techniques for gradient estimates in reinforcement learning. Journal of Machine Learning Research 5 (Nov), p. 1471–1530. Cited by: item ➊, §2.1. [47] W. W. Hager (1989) Updating the inverse of a matrix. SIAM review 31 (2), p. 221–239. Cited by: Lemma 1. [48] J. Hao, K. Ji, and M. Liu (2024) Bilevel coreset selection in continual learning: a new formulation and algorithm. Advances in Neural Information Processing Systems 36. Cited by: Appendix A, item ➊, item ➍, §1, §1, §1, §2.1, §2.1, §2.2, §2, §5, §5. [49] A. Hashemi, M. Ghasemi, H. Vikalo, and U. Topcu (2019) Submodular observation selection and information gathering for quadratic models. In International Conference on Machine Learning, p. 2653–2662. Cited by: §B.1, §4.1. [50] C. He, H. Ye, L. Shen, and T. Zhang (2020) Milenas: efficient neural architecture search via mixed-level reformulation. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 11993–12002. Cited by: Appendix A. [51] K. He, X. Zhang, S. Ren, and J. Sun (2016) Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, p. 770–778. Cited by: Appendix H, §5. [52] X. He, L. Chen, D. Cheng, V. Gupta, et al. (2022) Revisit last-iterate convergence of msgd under milder requirement on step size. Advances in Neural Information Processing Systems 35, p. 36559–36570. Cited by: Appendix L. [53] M. Hong, H. Wai, Z. Wang, and Z. Yang (2023) A two-timescale stochastic algorithm framework for bilevel optimization: complexity analysis and application to actor-critic. SIAM Journal on Optimization 33 (1), p. 147–180. Cited by: §2.1. [54] R. A. Horn and C. R. Johnson (2012) Matrix analysis. Cambridge university press. Cited by: Lemma 1. [55] X. Huang, Z. Liu, S. Liu, and K. Cheng (2024) Robust and efficient quantization-aware training via coreset selection. Transactions on Machine Learning Research. Cited by: §1. [56] K. Ji, J. Yang, and Y. Liang (2021) Bilevel optimization: convergence analysis and enhanced design. In International conference on machine learning, p. 4882–4892. Cited by: Appendix A, §2.1. [57] L. V. Jospin, H. Laga, F. Boussaid, W. Buntine, and M. Bennamoun (2022) Hands-on bayesian neural networks—a tutorial for deep learning users. IEEE Computational Intelligence Magazine 17 (2), p. 29–48. Cited by: §3.1. [58] S. M. Kakade (2001) A natural policy gradient. Advances in neural information processing systems 14. Cited by: §1. [59] G. Khashayar and G. R. Manuel (2019) Non-submodular function maximization subject to a matroid constraint, with applications. External Links: 1811.07863, Link Cited by: Appendix A. [60] K. Killamsetty, X. Zhao, F. Chen, and R. Iyer (2021) Retrieve: coreset selection for efficient and robust semi-supervised learning. Advances in neural information processing systems 34, p. 14488–14501. Cited by: §1. [61] J. Kirkpatrick, R. Pascanu, N. Rabinowitz, J. Veness, G. Desjardins, A. A. Rusu, K. Milan, J. Quan, T. Ramalho, A. Grabska-Barwinska, et al. (2017) Overcoming catastrophic forgetting in neural networks. Proceedings of the national academy of sciences 114 (13), p. 3521–3526. Cited by: Appendix A. [62] A. Krizhevsky G. Hinton et al. (2009) Learning multiple layers of features from tiny images. Cited by: §5. [63] L. Kumari, S. Wang, A. Das, T. Zhou, and J. Bilmes (2024) An end-to-end submodular framework for data-efficient in-context learning. In Findings of the Association for Computational Linguistics: NAACL 2024, p. 3293–3308. Cited by: §1. [64] L. Kumari, S. Wang, T. Zhou, N. Sarda, A. Rowe, and J. Bilmes (2024) BumbleBee: dynamic KV-cache streaming submodular summarization for infinite-context transformers. In First Conference on Language Modeling, Cited by: §1. [65] G. Lan (2020) First-order and stochastic optimization methods for machine learning. Vol. 1, Springer. Cited by: §1. [66] Z. Li and D. Hoiem (2017) Learning without forgetting. IEEE transactions on pattern analysis and machine intelligence 40 (12), p. 2935–2947. Cited by: Appendix A. [67] X. Lin, W. Wang, Y. Li, S. Yang, F. Feng, Y. Wei, and T. Chua (2024) Data-efficient fine-tuning for llm-based recommendation. In Proceedings of the 47th international ACM SIGIR conference on research and development in information retrieval, p. 365–374. Cited by: §1. [68] B. Liu, M. Ye, S. Wright, P. Stone, and Q. Liu (2022) Bome! bilevel optimization made easy: a simple first-order approach. Advances in neural information processing systems 35, p. 17248–17262. Cited by: Appendix A. [69] C. Liu, B. Zoph, M. Neumann, J. Shlens, W. Hua, L. Li, L. Fei-Fei, A. Yuille, J. Huang, and K. Murphy (2018) Progressive neural architecture search. In Proceedings of the European conference on computer vision (ECCV), p. 19–34. Cited by: Appendix A. [70] Y. Liu, W. Dong, L. Zhang, D. Gong, and Q. Shi (2019) Variational bayesian dropout with a hierarchical prior. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 7124–7133. Cited by: §3.1. [71] Z. Liu and Z. Zhou (2024) Revisiting the last-iterate convergence of stochastic gradient methods. In The Twelfth International Conference on Learning Representations, External Links: Link Cited by: Appendix L. [72] D. Lopez-Paz and M. Ranzato (2017) Gradient episodic memory for continual learning. Advances in neural information processing systems 30. Cited by: Appendix A, item ➋, §1. [73] J. Lorraine, P. Vicol, and D. Duvenaud (2020) Optimizing millions of hyperparameters by implicit differentiation. In International conference on artificial intelligence and statistics, p. 1540–1552. Cited by: §1. [74] Z. Lu and S. Mei (2024) First-order penalty methods for bilevel optimization. SIAM Journal on Optimization 34 (2), p. 1937–1969. Cited by: Appendix A. [75] D. Maclaurin, D. Duvenaud, and R. Adams (2015) Gradient-based hyperparameter optimization through reversible learning. In International conference on machine learning, p. 2113–2122. Cited by: Appendix A. [76] P. Manupriya, P. Jawanpuria, K. S. Gurumoorthy, S. Jagarlapudi, and B. Mishra (2024) Submodular framework for structured-sparse optimal transport. In Forty-first International Conference on Machine Learning, Cited by: §4.1. [77] M. McCloskey and N. J. Cohen (1989) Catastrophic interference in connectionist networks: the sequential learning problem. In Psychology of learning and motivation, Vol. 24, p. 109–165. Cited by: Appendix A, §1. [78] S. Mindermann, J. M. Brauner, M. T. Razzak, M. Sharma, A. Kirsch, W. Xu, B. Höltgen, A. N. Gomez, A. Morisot, S. Farquhar, et al. (2022) Prioritized training on points that are learnable, worth learning, and not yet learnt. In International Conference on Machine Learning, p. 15630–15649. Cited by: Appendix A, §1, §3.1, §3.1. [79] B. Mirzasoleiman, J. Bilmes, and J. Leskovec (2020) Coresets for data-efficient training of machine learning models. In International Conference on Machine Learning, p. 6950–6960. Cited by: §1. [80] B. Mirzasoleiman, A. Karbasi, R. Sarkar, and A. Krause (2013) Distributed submodular maximization: identifying representative elements in massive data. Advances in Neural Information Processing Systems 26. Cited by: §5. [81] K. P. Murphy (2012) Machine learning: a probabilistic perspective. The MIT Press. External Links: ISBN 0262018020 Cited by: §3.1. [82] B. K. Natarajan (1995) Sparse approximate solutions to linear systems. SIAM journal on computing 24 (2), p. 227–234. Cited by: §2.1. [83] Y. Nesterov (2013) Introductory lectures on convex optimization: a basic course. Vol. 87, Springer Science & Business Media. Cited by: §1. [84] Y. Nesterov and V. Spokoiny (2017) Random gradient-free minimization of convex functions. Foundations of Computational Mathematics 17 (2), p. 527–566. Cited by: §3.1. [85] C. V. Nguyen, Y. Li, T. D. Bui, and R. E. Turner (2017) Variational continual learning. arXiv preprint arXiv:1710.10628. Cited by: §5. [86] L. N. Nguyen and M. T. Thai (2022) Efficient algorithms for monotone non-submodular maximization with partition matroid constraint. In Proceedings of the Thirty-First International Joint Conference on Artificial Intelligence, IJCAI-22, L. D. Raedt (Ed.), p. 4807–4813. Note: Main Track Cited by: §4.1. [87] O. Ostapenko, P. Rodriguez, M. Caccia, and L. Charlin (2021) Continual learning via local module composition. Advances in Neural Information Processing Systems 34, p. 30298–30312. Cited by: Appendix A. [88] P. M. Pardalos and N. Kovoor (1990) An algorithm for a singly constrained class of quadratic programs subject to upper and lower bounds. Mathematical Programming 46, p. 321–328. Cited by: Appendix E. [89] F. Pedregosa (2016) Hyperparameter optimization with approximate gradient. In International conference on machine learning, p. 737–746. Cited by: Appendix A, §1. [90] T. Rainforth, A. Foster, D. R. Ivanova, and F. Bickford Smith (2024) Modern bayesian experimental design. Statistical Science 39 (1), p. 100–114. Cited by: §B.1, §4.1. [91] A. Rajeswaran, C. Finn, S. M. Kakade, and S. Levine (2019) Meta-learning with implicit gradients. Advances in neural information processing systems 32. Cited by: Appendix A. [92] S. Rebuffi, A. Kolesnikov, G. Sperl, and C. H. Lampert (2017) Icarl: incremental classifier and representation learning. In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, p. 2001–2010. Cited by: §5. [93] M. Riemer, I. Cases, R. Ajemian, M. Liu, I. Rish, Y. Tu, and G. Tesauro (2018) Learning to learn without forgetting by maximizing transfer and minimizing interference. arXiv preprint arXiv:1810.11910. Cited by: Appendix H. [94] H. Ritter, A. Botev, and D. Barber (2018) Online structured laplace approximations for overcoming catastrophic forgetting. Advances in Neural Information Processing Systems 31. Cited by: Appendix A. [95] J. Schreiber, J. Bilmes, and W. S. Noble (2020) Apricot: submodular selection for data summarization in python. Journal of Machine Learning Research 21 (161), p. 1–6. External Links: Link Cited by: §5. [96] O. Sener and S. Savarese (2017) Active learning for convolutional neural networks: a core-set approach. arXiv preprint arXiv:1708.00489. Cited by: §5. [97] A. Shaban, C. Cheng, N. Hatch, and B. Boots (2019) Truncated back-propagation for bilevel optimization. In The 22nd international conference on artificial intelligence and statistics, p. 1723–1732. Cited by: Appendix A. [98] O. Shamir (2017) An optimal algorithm for bandit and zero-order convex optimization with two-point feedback. Journal of Machine Learning Research 18 (52), p. 1–11. External Links: Link Cited by: §3.1. [99] D. Sow, K. Ji, and Y. Liang (2022) On the convergence theory for hessian-free bilevel algorithms. Advances in Neural Information Processing Systems 35, p. 4136–4149. Cited by: Appendix A. [100] S. Sujit, S. Nath, P. Braga, and S. Ebrahimi Kahou (2023) Prioritizing samples in reinforcement learning with reducible loss. Advances in Neural Information Processing Systems 36, p. 23237–23258. Cited by: §C.3. [101] R. S. Sutton, D. McAllester, S. Singh, and Y. Mansour (1999) Policy gradient methods for reinforcement learning with function approximation. Advances in neural information processing systems 12. Cited by: §C.1, item ➊, §1, §2.1. [102] T. Thiery and J. Ward (2022) Two-sided weak submodularity for matroid constrained optimization and regression. In Proceedings of Thirty Fifth Conference on Learning Theory, P. Loh and M. Raginsky (Eds.), Vol. 178, p. 3605–3634. Cited by: Appendix A, Appendix D, item ∙ , §1, §4.1. [103] S. Tofigh, M. Askarizadeh, M. O. Ahmad, M. Swamy, and K. K. Nguyen (2025) ACLI: a cnn pruning framework leveraging adjacent convolutional layer interdependence and γ-weakly submodularity. IEEE Transactions on Pattern Analysis and Machine Intelligence. Cited by: §4.1. [104] R. Tong, Y. Liu, J. Q. Shi, and D. Gong (2025) Coreset selection via reducible loss in continual learning. In The Thirteenth International Conference on Learning Representations, Cited by: Appendix A, Figure 5, Figure 5, 3rd item, 3rd item, §K.1, §K.1, §K.2, §K.2, §K.2, Appendix H, item ➋, §1, §2.1, §3.1, §5, §5, §5. [105] F. Wang, D. Zhou, H. Ye, and D. Zhan (2022) Foster: feature boosting and compression for class-incremental learning. In European conference on computer vision, p. 398–414. Cited by: Appendix A. [106] V. D. Wild, S. Ghalebikesabi, D. Sejdinovic, and J. Knoblauch (2023) A rigorous link between deep ensembles and (variational) bayesian methods. Advances in Neural Information Processing Systems 36, p. 39782–39811. Cited by: §3.1. [107] R. J. Williams (1992) Simple statistical gradient-following algorithms for connectionist reinforcement learning. Machine learning 8 (3), p. 229–256. Cited by: item ➊, §1, §2.1. [108] C. Wu, A. Rajeswaran, Y. Duan, V. Kumar, A. M. Bayen, S. Kakade, I. Mordatch, and P. Abbeel (2018) Variance reduction for policy gradient with action-dependent factorized baselines. In International Conference on Learning Representations, External Links: Link Cited by: item ➊, §2.1. [109] Q. Yan, D. Gong, Y. Liu, A. Van Den Hengel, and J. Q. Shi (2022) Learning bayesian sparse networks with full experience replay for continual learning. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 109–118. Cited by: Appendix A. [110] S. Yan, J. Xie, and X. He (2021) Der: dynamically expandable representation for class incremental learning. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, p. 3014–3023. Cited by: Appendix A. [111] J. Yoon, D. Madaan, E. Yang, and S. J. Hwang (2022) Online coreset selection for rehearsal-based continual learning. In International Conference on Learning Representations, External Links: Link Cited by: Appendix A, §1, §5. [112] F. Zenke, B. Poole, and S. Ganguli (2017) Continual learning through synaptic intelligence. In International conference on machine learning, p. 3987–3995. Cited by: Appendix A, §K.2, §5. [113] J. Zhang, Y. Qin, R. Pi, W. Zhang, R. Pan, and T. Zhang (2025) Tagcos: task-agnostic gradient clustered coreset selection for instruction tuning data. In Findings of the Association for Computational Linguistics: NAACL 2025, p. 4671–4686. Cited by: §1. [114] Y. Zhang, P. Khanduri, I. Tsaknakis, Y. Yao, M. Hong, and S. Liu (2024) An introduction to bilevel optimization: foundations and applications in signal processing and machine learning. IEEE Signal Processing Magazine 41 (1), p. 38–59. Cited by: Appendix A, item ➊. [115] B. Zhou, R. Jiang, and S. Shen (2024) Learning to solve bilevel programs with binary tender. In The Twelfth International Conference on Learning Representations, External Links: Link Cited by: §2.1. [116] X. Zhou, R. Pi, W. Zhang, Y. Lin, Z. Chen, and T. Zhang (2022) Probabilistic bilevel coreset selection. In International Conference on Machine Learning, p. 27287–27302. Cited by: Appendix A, Figure 5, Figure 5, 2nd item, 2nd item, §K.1, §K.1, §K.2, §C.1, item ➍, §1, §1, §1, §2.1, §2.1, §2.2, §2, §3, §5, §5, §5. Appendix Contents A Related Work .A B Proof of Theorem 2 and the Equivalence of the V-Optimal Criterion .B B.1 Proof of Theorem 2 .B B.2 Equivalence of the V-Optimal Criterion .B C Gradient Estimators for Φ⁡(π) ( π) .C.1 C.1 Policy Gradient Estimator .C.1 C.2 Two-point Gradient Estimator .C.2 C.3 Reducible Two-point Representation .C.3 D Proof of Theorem 3 .D E Projection Calculation .E F Effective Reservoir Rounding .F G Data Summarization on Cifar-100 Dataset .G H Coreset Selection for Existing Continual Learning Methods .H I Experiment with Adversarial Initialization .I J Data Corruption .J K Experimental Details .K K.1 Data Summarization Experiment Details .K.1 K.2 Continual Learning Experiment Details .K.2 L Limitations and Broader Impact .L Appendix A Related Work Bilevel Continuous Optimization. Bilevel continuous optimization (BCO) provides a natural framework for problems whose objective depends on the solution of another optimization problem, and has been widely used in hyperparameter optimization [37, 38, 97], meta-learning [36, 91], and neural architecture search [31, 69, 50]. Since its early formulation in mathematical programming [10], BCO has attracted extensive attention in both optimization and machine learning [114]. A major line of work focuses on hypergradient-based methods, which can be broadly categorized into approximate implicit differentiation (AID) methods [42, 44, 45, 56, 89] and iterative differentiation (ITD) methods [37, 56, 75]. AID methods compute the hypergradient by differentiating the optimality conditions of the lower-level problem, typically involving Hessian inverse-vector products or their approximations, while ITD methods differentiate through the trajectory of the lower-level optimization procedure. Another line of work adopts penalty-based reformulations, which transform the bilevel problem into a constrained or penalized single-level problem by replacing the lower-level solution requirement with its optimality conditions [21, 74, 99, 68]. Continual Learning. Continual learning aims to train a model on a sequence of tasks while preserving the knowledge acquired from previously seen data. A central challenge in this setting is catastrophic forgetting [39, 77], where the model quickly loses its performance on old tasks after being updated on new ones. Existing methods address this issue from different perspectives, including regularization-based methods [16, 61, 66, 94, 109, 112], parameter isolation based methods [87, 105, 110], and replay-based methods [8, 9, 12, 72, 104, 111, 116]. Among them, replay-based continual learning combats forgetting by maintaining a small memory buffer of representative samples from previous tasks and replaying them during future training. Since the memory budget is typically much smaller than the accumulated training data, the quality of the retained samples becomes critical to the effectiveness of rehearsal. This has motivated a series of coreset selection methods for continual learning [8, 9, 116, 48, 104]. For example, [8] extends coreset selection to deep neural networks by formulating it as a bilevel optimization problem. After that, [116] further introduces probabilistic masks to perform global subset selection. Next, [48] revisits bilevel coreset selection in continual learning and proposes a new 0-1 relaxation formulation. More recently, [104] design a greedy method according to the reducible loss proposed by [78]. Weakly Submodular Maximization. The notion of γ-weak submodularity was originally introduced by [24], who showed that the standard greedy algorithm achieves a (1−e−γ)(1-e^-γ)-approximation for maximizing a monotone γ-weakly submodular function under a cardinality constraint. Subsequently, [59] studied weakly submodular maximization beyond simple cardinality constraints and proved that the standard greedy algorithm achieves a 0.4γ2rγ+1 0.4γ^2 rγ+1-approximation for maximizing a monotone γ-weakly submodular function under a matroid constraint, where r denotes the rank of the matroid. Unfortunately, this approximation guarantee explicitly depends on r and therefore deteriorates as the matroid rank increases. Note that the partition constraints considered in this paper are a special case of matroid constraints. Later, [102] introduced the upper submodularity ratio β and developed a distorted local search algorithm for monotone (γ,β)(γ,β)-weakly submodular maximization. This distorted local search typically requires repeated swap operations and will incur substantial computational overhead, making it difficult to apply to large-scale selection problems. Appendix B Proof of Theorem 2 and the Equivalence of the V-Optimal Criterion B.1 Proof of Theorem 2 In this subsection, we show the detailed verification of Theorem 2. Before that, we recall the bilevel data selection problem for Bayesian linear regression. Specifically, the data is generated from a linear model ≜X+ϵ y X θ+ ε, ϵ ε, where ≜[y1,…,yn]⊤∈ℝn y [y_1,…,y_n] ^n, X≜[x1,…,xn]⊤∈ℝn×pX [x_1,…,x_n] ^n× p, the noise vector ϵ ε follows a Gaussian distribution ⁡(,σ2In)N( 0,σ^2I_n) and the model parameter θ is equipped with a Gaussian prior, i.e., ∼⁡(,Λ) θ ( 0, ). In this case, the lower-level optimization of bilevel data selection problem can be written as ∗()≜arg⁡min⁡∑(xi,yi)∈(yi−xi⊤)2+σ2⊤Λ−1 θ^*(S) _ θ _(x_i,y_i) (y_i-x_i θ)^2+σ^2 θ ^-1 θ, and our goal is to select a small subset ⊆≜x1,…,xnS \x_1,…,x_n\ such that the resulting posterior estimator ∗() θ^*(S) has a small prediction error over the full dataset. More importantly, from Bayesian statistical theory [14, 49, 90], we can show that ∗() θ^*(S) has the following closed form: ∗()≜(XTX+σ2Λ−1)−1XT, θ^*(S) (X_S^TX_S+σ^2 ^-1)^-1X_S^T y_S, where X_S and y_S denotes the matrix X and the vector y restricted to rows induced by S, respectively. Since the proof of Theorem 2 will repeatedly involve matrix inverses, we begin by recalling the classical Sherman-Morrisson-Woodbury formula, that is to say, Lemma 1 (Sherman-Morrisson-Woodbury formula [54, 47]). Let A, U, C and V be matrices of conformable sizes. Then, (A+UCV)−1≜A−1−A−1U(C−1+VA−1U)−1VA−1(A+UCV)^-1 A^-1-A^-1U(C^-1+VA^-1U)^-1VA^-1. With all above concepts and formula, we next prove Theorem 2. Proof. We first simplify the Bayesian V-optimal criterion ℒV(∗())≜12nϵ,[‖X⁡(−∗())‖22]L_V( θ^*(S)) 12nE_ ε, θ[\|X ( θ- θ^*(S) )\|_2^2], namely, ℒV(∗())≜12nϵ,[‖X⁡(∗()−)‖22] _V( θ^*(S)) 12nE_ ε, θ[\|X ( θ^*(S)- θ )\|_2^2] =12nϵ,[‖X⁡((XTX+σ2Λ−1)−1XT−)‖22] = 12nE_ ε, θ[\|X ((X_S^TX_S+σ^2 ^-1)^-1X_S^T y_S- θ )\|_2^2] =12nϵ,[‖X⁡((XTX+σ2Λ−1)−1XT(X+ϵ)−)‖22] = 12nE_ ε, θ[\|X ((X_S^TX_S+σ^2 ^-1)^-1X_S^T (X_S θ+ ε_S )- θ )\|_2^2] =12nϵ,[∥X((XTX+σ2Λ−1)−1XTX−+X((XTX+σ2Λ−1)−1XTϵ)∥22] = 12nE_ ε, θ[\|X ((X_S^TX_S+σ^2 ^-1)^-1X_S^TX_S θ- θ+X ((X_S^TX_S+σ^2 ^-1)^-1X_S^T ε_S )\|_2^2] =12nϵ,[‖X⁡((XTX+σ2Λ−1)−1XTϵ−σ2(XTX+σ2Λ−1)−1Λ−1)‖22] = 12nE_ ε, θ[\|X ((X_S^TX_S+σ^2 ^-1)^-1X_S^T ε_S-σ^2(X_S^TX_S+σ^2 ^-1)^-1 ^-1 θ )\|_2^2] =12nϵ,[∥X(XTX+σ2Λ−1)−1XTϵ∥22]+12nϵ,[∥σ2⋅X(XTX+σ2Λ−1)−1Λ−1)∥22] = 12nE_ ε, θ[\|X(X_S^TX_S+σ^2 ^-1)^-1X_S^T ε_S\|_2^2]+ 12nE_ ε, θ[\|σ^2· X(X_S^TX_S+σ^2 ^-1)^-1 ^-1 θ )\|_2^2] =12nTr[X(XTX+σ2Λ−1)−1XTϵ[ϵϵT]X(XTX+σ2Λ−1)−1XT] = 12nTr[X(X_S^TX_S+σ^2 ^-1)^-1X_S^TE_ ε [ ε_S ε_S^T ]X_S(X_S^TX_S+σ^2 ^-1)^-1X^T] +12nTr[σ2⋅X(XTX+σ2Λ−1)−1Λ−1[T]Λ−1(XTX+σ2Λ−1)−1XT⋅σ2] + 12nTr[σ^2· X(X_S^TX_S+σ^2 ^-1)^-1 ^-1E_ θ [ θ θ^T ] ^-1(X_S^TX_S+σ^2 ^-1)^-1X^T·σ^2] =σ22nTr[X(XTX+σ2Λ−1)−1XTX(XTX+σ2Λ−1)−1XT] = σ^22nTr[X_S(X_S^TX_S+σ^2 ^-1)^-1X^TX(X_S^TX_S+σ^2 ^-1)^-1X_S^T] +σ22nTr[σ2⋅X(XTX+σ2Λ−1)−1Λ−1(XTX+σ2Λ−1)−1XT] + σ^22nTr[σ^2· X(X_S^TX_S+σ^2 ^-1)^-1 ^-1(X_S^TX_S+σ^2 ^-1)^-1X^T] =σ22nTr[(XTX+σ2Λ−1)−1XTX(XTX+σ2Λ−1)−1XTX] = σ^22nTr[(X_S^TX_S+σ^2 ^-1)^-1X^TX(X_S^TX_S+σ^2 ^-1)^-1X_S^TX_S] +σ22nTr[(XTX+σ2Λ−1)−1XTX(XTX+σ2Λ−1)−1(σ2⋅Λ−1)] + σ^22nTr[(X_S^TX_S+σ^2 ^-1)^-1X^TX(X_S^TX_S+σ^2 ^-1)^-1 (σ^2· ^-1 )] =σ22nTr[(XTX+σ2Λ−1)−1XTX]=σ22nTr[X(XTX+σ2Λ−1)−1XT], = σ^22nTr[(X_S^TX_S+σ^2 ^-1)^-1X^TX]= σ^22nTr[X(X_S^TX_S+σ^2 ^-1)^-1X^T], where ϵ[ϵϵT]=σ2I||E_ ε [ ε_S ε_S^T ]=σ^2I_|S| and [T]=ΛE_ θ [ θ θ^T ]= . With this equality, we can rewrite the utility function UV()U_V(S) as UV() U_V(S) ≜ℒV()−ℒV(∗()) _V( 0)-L_V( θ^*(S)) =σ22n(Tr⁡[σ−2⋅XΛXT]−Tr⁡[X(XTX+σ2Λ−1)−1XT]). = σ^22n (Tr[σ^-2· X X^T]-Tr[X(X_S^TX_S+σ^2 ^-1)^-1X^T] ). For simplicity, in the remainder of this proof, for any subset S, we set ℳ≜XTX+σ2Λ−1M_S X_S^TX_S+σ^2 ^-1. Next, we show this utility function UV()U_V(S) satisfies the (γ,β)(γ,β)-weak submodularity. At first, for any subset ⊆ℬA and ℛ≜ℬ∖R , we have that UV(ℬ)−UV() U_V(B)-U_V(A) =σ22n(Tr[Xℳ−1]XT]−Tr[Xℳ−1ℬ]XT]) = σ^22n (Tr[XM^-1_A]X^T]-Tr[XM^-1_B]X^T] ) =σ22n(Tr[Xℳ−1]XT]−Tr[X(ℳ+XℛTXℛ)−1XT]) = σ^22n (Tr[XM^-1_A]X^T]-Tr[X (M_A+X_R^TX_R )^-1X^T] ) =σ22n(Tr[Xℳ−1]XT]−Tr[X(ℳ−1−ℳ−1XℛT(I+Xℛℳ−1XℛT)−1Xℛℳ−1)XT]) = σ^22n (Tr[XM^-1_A]X^T]-Tr[X (M_A^-1-M_A^-1X_R^T (I+X_RM_A^-1X_R^T )^-1X_RM_A^-1 )X^T] ) =σ22n(Tr⁡[Xℳ−1XℛT(I+Xℛℳ−1XℛT)−1Xℛℳ−1XT]) = σ^22n (Tr[XM_A^-1X_R^T (I+X_RM_A^-1X_R^T )^-1X_RM_A^-1X^T] ) =σ22n(Tr⁡[(I+Xℛℳ−1XℛT)−1Xℛℳ−1XTXℳ−1XℛT]) = σ^22n (Tr[ (I+X_RM_A^-1X_R^T )^-1X_RM_A^-1X^TXM_A^-1X_R^T] ) Note that λmin(I+Xℛℳ−1XℛT)≥1 _min(I+X_RM_A^-1X_R^T)≥ 1, so we have that UV(ℬ)−UV()≤σ22nTr[Xℛℳ−1XTXℳ−1XℛT].U_V(B)-U_V(A)≤ σ^22nTr[X_RM_A^-1X^TXM_A^-1X_R^T]. Furthermore, we also have that (Note that xe∈ℝp×1x_e ^p× 1) UV(∪e)−UV() U_V(A∪\e\)-U_V(A) =σ22n(Tr⁡[(1+xeTℳ−1xe)−1xeTℳ−1XTXℳ−1xe]) = σ^22n (Tr[ (1+x_e^TM_A^-1x_e )^-1x^T_eM_A^-1X^TXM_A^-1x_e] ) =σ22n⋅Tr⁡[xeTℳ−1XTXℳ−1xe]1+xeTℳ−1xe = σ^22n· Tr[x_e^TM_A^-1X^TXM_A^-1x_e]1+x_e^TM_A^-1x_e ≥σ22n⋅Tr⁡[xeTℳ−1XTXℳ−1xe]1+s2σ2⋅λmax(Λ), ≥ σ^22n· Tr[x_e^TM_A^-1X^TXM_A^-1x_e]1+ s^2σ^2· _max( ), where the final inequality follows from ℳ=XTX+σ2Λ−1M_A=X_S^TX_S+σ^2 ^-1 such that xeTℳ−1xe≤xeT(σ2Λ−1)−1xe≤xeTΛxeσ2≤s2σλmax(Λ)x_e^TM_A^-1x_e≤ x_e^T (σ^2 ^-1 )^-1x_e≤ x_e^T x_eσ^2≤ s^2σ _max( ). Here, s≜maxi∈[n]⁡‖xi‖2s _i∈[n]\|x_i\|_2 and λmax(Λ) _ ( ) is the largest eigenvalue of the matrix Λ . As a result, we have that ∑e∈ℬ∖(UV(∪e)−UV()) _e (U_V(A∪\e\)-U_V(A) ) ≥σ22n⋅∑e∈ℬ∖Tr⁡[xeTℳ−1XTXℳ−1xe]1+s2σ2⋅λmax(Λ) ≥ σ^22n· _e Tr[x_e^TM_A^-1X^TXM_A^-1x_e]1+ s^2σ^2· _max( ) =σ22n⋅Tr⁡[Xℛℳ−1XTXℳ−1XℛT]1+s2σ2⋅λmax(Λ)≥UV(ℬ)−UV()1+s2σ2⋅λmax(Λ), = σ^22n· Tr[X_RM_A^-1X^TXM_A^-1X_R^T]1+ s^2σ^2· _max( )≥ U_V(B)-U_V(A)1+ s^2σ^2· _max( ), where ℛ≜ℬ∖R . Secondly, for any subset ⊆ℬA and ℛ≜ℬ∖R , we have that UV(ℬ)−UV() U_V(B)-U_V(A) =σ22n(Tr[Xℳ−1]XT]−Tr[Xℳ−1ℬ]XT]) = σ^22n (Tr[XM^-1_A]X^T]-Tr[XM^-1_B]X^T] ) =σ22n(Tr[X(ℳℬ−XℛTXℛ)−1XT]−Tr[Xℳ−1ℬ]XT]) = σ^22n (Tr[X (M_B-X_R^TX_R )^-1X^T]-Tr[XM^-1_B]X^T] ) =σ22n(Tr[X(ℳℬ−1+ℳℬ−1XℛT(I−Xℛℳℬ−1XℛT)−1Xℛℳℬ−1)XT−Tr[Xℳ−1ℬ]XT]]) = σ^22n (Tr[X (M_B^-1+M_B^-1X_R^T (I-X_RM_B^-1X_R^T )^-1X_RM_B^-1 )X^T-Tr[XM^-1_B]X^T]] ) =σ22n(Tr⁡[Xℳℬ−1XℛT(I−Xℛℳℬ−1XℛT)−1Xℛℳℬ−1XT]) = σ^22n (Tr[XM_B^-1X_R^T (I-X_RM_B^-1X_R^T )^-1X_RM_B^-1X^T] ) =σ22n(Tr⁡[(I−Xℛℳℬ−1XℛT)−1Xℛℳℬ−1XTXℳℬ−1XℛT]) = σ^22n (Tr[ (I-X_RM_B^-1X_R^T )^-1X_RM_B^-1X^TXM_B^-1X_R^T] ) Note that λmax(I−Xℛℳℬ−1XℛT)≤1 _max(I-X_RM_B^-1X_R^T)≤ 1, so we have that UV(ℬ)−UV()≥σ22nTr[Xℛℳℬ−1XTXℳℬ−1XℛT].U_V(B)-U_V(A)≥ σ^22nTr[X_RM_B^-1X^TXM_B^-1X_R^T]. Furthermore, we also have that UV(ℬ)−UV(ℬ∖e)=σ22n(Tr⁡[(1−xeTℳℬ−1xe)−1xeTℳℬ−1XTXℳℬ−1xe]). U_V(B)-U_V(B \e\)= σ^22n (Tr[ (1-x_e^TM_B^-1x_e )^-1x_e^TM_B^-1X^TXM_B^-1x_e] ). From Sherman-Morrisson-Woodbury formula, we also can show that (1−xeTℳℬ−1xe)−1=1+xeT(ℳℬ−xexeT)−1xe=1+xeTℳℬ∖e−1xe≤1+s2σλmax(Λ), (1-x_e^TM_B^-1x_e )^-1=1+x_e^T(M_B-x_ex_e^T)^-1x_e=1+x_e^TM_B \e\^-1x_e≤ 1+ s^2σ _max( ), where ℳℬ∖e≜Xℬ∖eTXℬ∖e+σ2Λ−1M_B \e\ X_B \e\^TX_B \e\+σ^2 ^-1 such that xeTℳℬ∖e−1xe≤xeT(σ2Λ−1)−1xe≤xeTΛxeσ2≤s2σλmax(Λ)x_e^TM_B \e\^-1x_e≤ x_e^T (σ^2 ^-1 )^-1x_e≤ x_e^T x_eσ^2≤ s^2σ _max( ). As a result, we have that UV(ℬ)−UV(ℬ∖e)≤σ22n(1+s2σλmax(Λ))⋅Tr⁡[xeTℳℬ−1XTXℳℬ−1xe].U_V(B)-U_V(B \e\)≤ σ^22n (1+ s^2σ _max( ) )·Tr[x_e^TM_B^-1X^TXM_B^-1x_e]. So we can have that ∑e∈ℬ∖(UV(ℬ)−UV(ℬ∖e)) _e (U_V(B)-U_V(B \e\) ) ≤σ22n(1+s2σλmax(Λ))⋅∑e∈ℬ∖Tr⁡[xeTℳℬ−1XTXℳℬ−1xe] ≤ σ^22n (1+ s^2σ _max( ) )· _e Tr[x_e^TM_B^-1X^TXM_B^-1x_e] =σ22n(1+s2σλmax(Λ))⋅Tr⁡[Xℛℳℬ−1XTXℳℬ−1XℛT] = σ^22n (1+ s^2σ _max( ) )·Tr[X_RM_B^-1X^TXM_B^-1X_R^T] ≤(1+s2σλmax(Λ))⋅(UV(ℬ)−UV()). ≤ (1+ s^2σ _max( ) )· (U_V(B)-U_V(A) ). In summary, The utility function UV()≜ℒV()−ℒV(∗())U_V(S) _V( 0)-L_V( θ^*(S)) be the induced by the Bayesian V-optimal criterion ℒV(∗())L_V( θ^*(S)) is monotone and (1/c,c)(1/c,c)-weakly submodular, where c≜1+s2σ2⋅λmax(Λ)c 1+ s^2σ^2· _max( ), s≜maxi∈[n]⁡‖xi‖2s _i∈[n]\|x_i\|_2 and λmax(Λ) _ ( ) is the largest eigenvalue of the matrix Λ . ∎ B.2 Equivalence of the V-Optimal Criterion In this subsection, we show the equivalence between the total prediction error ℒ⁡(∗())≜12nϵ,[‖−X∗()‖22]L( θ^*(S)) 12nE_ ε, θ [\| y-X θ^*(S)\|_2^2 ] and the Bayesian V-optimal criterion ℒV(∗())≜12nϵ,[‖X⁡(−∗())‖22]L_V( θ^*(S)) 12nE_ ε, θ[\|X ( θ- θ^*(S) )\|_2^2]. Specifically, when the number of data points n is sufficiently large, these two objectives differ only by an additive constant σ2/2σ^2/2. Consequently, the weak submodularity of ℒV(∗())L_V( θ^*(S)) established above, to some extent, also implies the weak submodularity of the total prediction error ℒ(∗())L( θ^*(S)). Lemma 2. Under the aforementioned assumptions of Bayesian linear regression, we can show that limn→∞ℒ⁡(∗())−ℒV(∗())=σ22. _n→∞L( θ^*(S))-L_V( θ^*(S))= σ^22. Proof. Fixed θ, we have, ℒ(∗()) ( θ^*(S)) = = 12nϵ[‖−X∗()‖22] 12nE_ ε [\| y-X θ^*(S)\|_2^2 ] = = 12nϵ[‖X−X∗()‖22]−1nϵ[ϵ⊤(X∗()−X)]+12nϵ[‖ϵ‖22] 12nE_ ε [\|X θ-X θ^*(S)\|^2_2 ]- 1nE_ ε [ε (X θ^*(S)-X θ) ]+ 12nE_ ε [\| ε\|_2^2 ] = = ℒV(∗())−1nϵ[ϵ⊤(X∗()−X)]+σ22 _V( θ^*(S))- 1nE_ ε [ε (X θ^*(S)-X θ) ]+ σ^22 = = ℒV(∗())−1nϵ[∑i∈Sϵi(xiT∗()−xiT)]−1nϵ[∑i∈[n]∖Sϵi(xiT∗()−xiT)]+σ22 _V( θ^*(S))- 1nE_ ε [ _i∈ S _i(x_i^T θ^*(S)-x_i^T θ) ]- 1nE_ ε [ _i∈[n] S _i(x_i^T θ^*(S)-x_i^T θ) ]+ σ^22 When n→∞n→∞, we have limn→∞1nϵ[∑i∈Sϵi(xiT∗()−xiT)]=0 _n→∞ 1nE_ ε [ _i∈ S _i(x_i^T θ^*(S)-x_i^T θ) ]=0 since S is a finite set. Since ∗() θ^*(S) only depends on ϵ ε_S, we have ϵ[∑i∈[n]∖Sϵi(xiT∗()−xiT)]=∑i∈[n]∖Sϵ[ϵi]ϵ[xiT∗()−xiT]=0.E_ ε [ _i∈[n] S _i(x_i^T θ^*(S)-x_i^T θ) ]= _i∈[n] SE_ ε [ _i ]E_ ε [x_i^T θ^*(S)-x_i^T θ ]=0. As a consequence, as limn→∞ℒ⁡(∗())−ℒV(∗())=σ22 _n→∞L( θ^*(S))-L_V( θ^*(S))= σ^22. ∎ Appendix C Gradient Estimators for Φ⁡() ( π) In this section, we will show two different gradient estimators for our proposed probabilistic relaxation Φ⁡() ( π), namely policy gradient estimator and two-point estimator in Eq.(8). C.1 Policy Gradient Estimator At first, from Eq.(7), we know that Φ⁡() ( π) ≜di,k∼i,∀i,k[ℒ(∗(∪i=1h∪k=1Kidi,k))] _d_i,k π_i,∀ i,k [L ( θ^* ( _i=1^h _k=1^K_id_i,k ) ) ] =∑di,k∈i∪∅,∀i,k(ℒ(∗(∪i=1h∪k=1Kidi,k))⋅∏i=1h∏k=1Kip(di,k∣i)), = _d_i,k _i∪\ \,∀ i,k (L ( θ^* ( _i=1^h _k=1^K_id_i,k ) )· _i=1^h _k=1^K_ip(d_i,k π_i) ), where the symbol di,k∼id_i,k π_i means that we sample di,kd_i,k from i∪∅D_i∪\ \ according to the categorical distribution i π_i. Then, like the policy gradient methods [116, 101, 4], we can show that ∂Φ⁡()∂πi,j=∑di,k∈i∪∅,∀i∈[h],k∈[Ki](ℒ(∗(∪i=1h∪k=1Kidi,k))∇πi,j(∏i=1h∏k=1Kip(di,k∣i))) ∂ ( π)∂ _i,j= _d_i,k _i∪\ \,∀ i∈[h],k∈[K_i] (L ( θ^* ( _i=1^h _k=1^K_id_i,k ) ) _ _i,j ( _i=1^h _k=1^K_ip(d_i,k π_i) ) ) =∑di,k∈i∪∅,∀i,k(ℒ(∗(∪i=1h∪k=1Kidi,k))∇πi,jlog(∏i=1h∏k=1Kip(di,k∣i))(∏i=1h∏k=1Kip(di,k∣i))) = _d_i,k _i∪\ \,∀ i,k (L ( θ^* ( _i=1^h _k=1^K_id_i,k ) ) _ _i,j ( _i=1^h _k=1^K_ip(d_i,k π_i) ) ( _i=1^h _k=1^K_ip(d_i,k π_i) ) ) =di,k∼i,∀i∈[h],k∈[Ki][ℒ(∗(∪i=1h∪k=1Kidi,k))∇πi,jlog(∏i=1h∏k=1Kip(di,k∣i))] =E_d_i,k π_i,∀ i∈[h],k∈[K_i] [L ( θ^* ( _i=1^h _k=1^K_id_i,k ) ) _ _i,j ( _i=1^h _k=1^K_ip(d_i,k π_i) ) ] =di,k∼i,∀i∈[h],k∈[Ki][ℒ(∗(∪i=1h∪k=1Kidi,k))∑k=1Ki([di,k=(xj(i),yj(i))]πi,j−[di,k=∅]1−∑k=1Kiπi,k)], =E_d_i,k π_i,∀ i∈[h],k∈[K_i] [L ( θ^* ( _i=1^h _k=1^K_id_i,k ) ) _k=1^K_i ( I[d_i,k=(x_j^(i),y_j^(i))] _i,j- I[d_i,k= ]1- _k=1^K_i _i,k ) ], where the second equality follows from dlog⁡(x)dx=1/x d (x)dx=1/x and ⁡[⋅]I[·] represents the indictor function. Note that, like Eq. (5), the policy gradient estimator of our proposed probabilistic relaxation still involves an unbounded score factor ([di,k=(xj(i),yj(i))]πi,j−[di,k=∅]1−∑k=1Kiπi,k) ( I[d_i,k=(x_j^(i),y_j^(i))] _i,j- I[d_i,k= ]1- _k=1^K_i _i,k ). Specifically, when πi,k→0 _i,k→ 0 or ∑k=1Kiπi,k→1 _k=1^K_i _i,k→ 1, the factor ([di,k=(xj(i),yj(i))]πi,j−[di,k=∅]1−∑k=1Kiπi,k) ( I[d_i,k=(x_j^(i),y_j^(i))] _i,j- I[d_i,k= ]1- _k=1^K_i _i,k ) will also become unbounded. C.2 Two-point Gradient Estimator In this subsection, we verify Theorem 1. Proof. At first, we can show that ∂Φ⁡()∂πi,j=∑di^,k^∈i^∪∅,∀i^∈[h],k^∈[Ki](ℒ(∗(∪i^=1h∪k^=1Kidi^,k^))∇πi,j(∏i^=1h∏k^=1Kip(di^,k^∣i^))) ∂ ( π)∂ _i,j= _d_ i, k _ i∪\ \,∀ i∈[h], k∈[K_i] (L ( θ^* ( _ i=1^h _ k=1^K_id_ i, k ) ) _ _i,j ( _ i=1^h _ k=1^K_ip(d_ i, k π_ i) ) ) (12) =∑di^,k^∈i^∪∅,∀i^,k^(ℒ(∗(∪i^=1h∪k^=1Kidi^,k^))∑k=1Ki(∇πi,jp(di,k∣i)⋅∏(i^,k^)≠(i,k)p(di^,k^∣i^))) = _d_ i, k _ i∪\ \,∀ i, k (L ( θ^* ( _ i=1^h _ k=1^K_id_ i, k ) ) _k=1^K_i ( _ _i,jp(d_i,k π_i)· _( i, k)≠(i,k)p(d_ i, k π_ i) ) ) =∑k=1Ki(∑di^,k^∈i^∪∅,∀i^,k^(ℒ(∗(∪i^=1h∪k^=1Kidi^,k^))∇πi,jp(di,k∣i)⋅∏(i^,k^)≠(i,k)p(di^,k^∣i^))) = _k=1^K_i ( _d_ i, k _ i∪\ \,∀ i, k (L ( θ^* ( _ i=1^h _ k=1^K_id_ i, k ) ) _ _i,jp(d_i,k π_i)· _( i, k)≠(i,k)p(d_ i, k π_ i) ) ) =Ki⋅(∑di^,k^∈i^∪∅,∀i^,k^(ℒ(∗(∪i^=1h∪k^=1Kidi^,k^))∇πi,jp(di,1∣i)⋅∏(i^,k^)≠(i,1)p(di^,k^∣i^))) =K_i· ( _d_ i, k _ i∪\ \,∀ i, k (L ( θ^* ( _ i=1^h _ k=1^K_id_ i, k ) ) _ _i,jp(d_i,1 π_i)· _( i, k)≠(i,1)p(d_ i, k π_ i) ) ) =Ki(∑di^,k^∈i^∪∅,(i^,k^)≠(i,1)(∑di,1∈i∪∅ℒ(∗(∪i^,k^di^,k^))∇πi,jp(di,1∣i))∏(i^,k^)≠(i,1)p(di^,k^∣i^)), =K_i ( _d_ i, k _ i∪\ \,( i, k)≠(i,1) ( _d_i,1 _i∪\ \L ( θ^* ( _ i, kd_ i, k ) ) _ _i,jp(d_i,1 π_i) ) _( i, k)≠(i,1)p(d_ i, k π_ i) ), where the fourth equality comes from the independence between di,1,di,2,…,di,Ki\d_i,1,d_i,2,…,d_i,K_i\. When di,1=(xj(i),yj(i))d_i,1=(x_j^(i),y_j^(i)) or di,1=∅d_i,1= , we can know that ∇πi,jp(di,1∣i)=1 _ _i,jp(d_i,1 π_i)=1 or ∇πi,jp(di,1∣i)=−1 _ _i,jp(d_i,1 π_i)=-1, respectively. As for other cases, we have ∇πi,jp(di,1∣i)=0 _ _i,jp(d_i,1 π_i)=0. Thus, if we define −(i,1)≜⋃(i′,k′)≠(i,1)di′,k′S_-(i,1) _(i ,k )≠(i,1)d_i ,k as the subset obtained by excluding the first draw di,1d_i,1 from source iD_i, we have the following equality: ∑di,1∈i∪∅ℒ(∗(∪i^,k^di,k))∇πi,jp(di,1∣i)=ℒ(∗(−(i,1)∪(x(i)j,y(i)j)))−ℒ(∗(−(i,1))). _d_i,1 _i∪\ \L ( θ^* ( _ i, kd_i,k ) ) _ _i,jp(d_i,1 π_i)=L ( θ^* (S_-(i,1)∪(x^(i)_j,y^(i)_j) ) )-L ( θ^* (S_-(i,1) ) ). As a result, we have that ∂Φ⁡()∂πi,j≜Ki⋅⁡[ℒ⁡(∗(−(i,1)∪(xj(i),yj(i))))−ℒ⁡(∗(−(i,1)))], ∂ ( π)∂ _i,j K_i·E [ ( θ^* (S_-(i,1)∪(x^(i)_j,y^(i)_j) ) )-L ( θ^* (S_-(i,1) ) ) ], (13) where the expectation is taken over all random draws except di,1d_i,1. ∎ C.3 Reducible Two-point Representation At first, we review the Bayes rule: for any random event A, B and C, we have p⁡(A∣B,C)≜p⁡(B∣A,C)⋅p⁡(A∣C)p⁡(B∣C).p(A B,C) p(B A,C)·p(A C)p(B C). As a result, in Eq.(9), if we set A=A= y, B=yB=y and C=,x,C=\ x,x,S\, we have that p⁡(|;∪(x,y)) ( y| x;S∪(x,y)) =p⁡(y|x;∪(,))⋅p⁡(|,x,)p⁡(y|x,,) = p(y|x;S∪( x, y))·p( y| x,x,S)p(y|x, x,S) =p⁡(y|x;∪)⋅p⁡(|,)p⁡(y|x,)=p⁡(y|x;)⋅p⁡(|,)p⁡(y|x,), = p(y|x;S )·p( y| x,S)p(y|x,S)= p(y|x;D)·p( y| x,S)p(y|x,S), where the second equality comes from the the conditional independence p⁡(yi|xi,xj;)=p⁡(yi|xi;)p(y_i|x_i,x_j;S)=p(y_i|x_i;S)(Please refer to Eq.(2) in [100]). As a result, we have that log⁡p⁡(|;∪(x,y))−log⁡p⁡(|,)=log⁡p⁡(y|x;)−log⁡p⁡(y|x,). ( y| x;S∪(x,y))- ( y| x,S)= (y|x;D)- (y|x,S). Furthermore, in Bayesian probabilistic framework, the loss function is commonly specified as the negative log-likelihood of the predictive distribution, namely, ℓ⁡(x,y,∗())≜−log⁡p⁡(y|x;) (x,y; θ^*(S)) - (y|x;S), then we can have that, throughout the independence between samples: ℒ⁡(∗())=∑i=1nℓ⁡(xi,yi,∗()) ( θ^* (S ) )= _i=1^n (x_i,y_i; θ^*(S)) =−∑i=1nlogp(yi|xi;)=−log(∏i=1np(yi|xi;))=−logp(|;). =- _i=1^n (y_i|x_i;S)=- ( _i=1^np(y_i|x_i;S) )=- ( y| x;S). As a result, we have ℒ⁡(∗(∪(x,y)))−ℒ⁡(∗()) ( θ^* (S∪(x,y) ) )-L ( θ^* (S ) ) =log⁡p⁡(|;)−log⁡p⁡(|;∪(x,y)) = ( y| x;S)- ( y| x;S∪(x,y)) =log⁡p⁡(y|x;)−log⁡p⁡(y|x,)=ℓ⁡(x,y,∗())−ℓ⁡(x,y,∗()). = (y|x;D)- (y|x,S)= (x,y; θ^*(D) )- (x,y; θ^*(S) ). Appendix D Proof of Theorem 3 In this section, we prove the core theorem of this paper, namely, Theorem 3. Before that, we prove a key first-order inequality regarding the upper-level utility function U⁡(∗())≜ℒ⁡(init)−ℒ⁡(∗())U( θ^*(S)) ( θ_init)-L( θ^*(S)). Next, we define some commonly used notations. At first, we set the corresponding probabilistic relaxation of he upper-level utility function U(∗())U( θ^*(S)) as Ψ⁡() ( π), that is, Ψ⁡() ( π) ≜∑di,k∈i∪∅,∀i,k(U(∗(∪i=1h∪k=1Kidi,k))⋅∏i=1h∏k=1Kip(di,k∣i)) _d_i,k _i∪\ \,∀ i,k (U ( θ^* ( _i=1^h _k=1^K_id_i,k ) )· _i=1^h _k=1^K_ip(d_i,k π_i) ) =∑di,k∈i∪∅,∀i,k((ℒ(init)−ℒ(∗(∪i=1h∪k=1Kidi,k)))⋅∏i=1h∏k=1Kip(di,k∣i)) = _d_i,k _i∪\ \,∀ i,k ( (L( θ_init)-L ( θ^* ( _i=1^h _k=1^K_id_i,k ) ) )· _i=1^h _k=1^K_ip(d_i,k π_i) ) =ℒ⁡(init)−Φ⁡(), =L( θ_init)- ( π), where Φ⁡() ( π) is our proposed probabilistic relaxation in Eq.(7). Furthermore, for any feasible data subset S satisfying |∩i|=Ki|S _i|=K_i for all i∈[h]i∈[h], we define its corresponding uniform prior as ≜∑i=1h1Ki∩i π_S _i=1^h 1K_i1_S _i. In other words, for each source iD_i, π_S places mass 1/Ki1/K_i on every candidate in ∩iS _i and mass 00 on the remaining candidates. Note that 1_A denotes the indicator vector of subset A, whose entry is 11 if data candidates is included in A and 00 otherwise. With the above preparations in place, we then show that Theorem 4. When the utility function U⁡(∗())≜ℒ⁡(init)−ℒ⁡(∗())U( θ^*(S)) ( θ_init)-L( θ^*(S)) induced by Problem Eq. (6) is (γ,β)(γ,β)-weakly monotone submodular, for any feasible subset ℬB, namely, |ℬ∩i|≤Ki,∀i∈[h]|B _i|≤ K_i,∀ i∈[h], the following inequality holds: γ2⋅U⁡(∗(ℬ))−(β+β⁡(1−γ)+γ2)⋅Ψ⁡()≤⟨ℬ−,∇Ψ()⟩.γ^2· U( θ^*(B))-(β+β(1-γ)+γ^2)· ( π)≤ π_B- π,∇ ( π) . Proof. At first, we investigate ⟨ℬ,∇Ψ()⟩ π_B,∇ ( π) . To this end, for each source iD_i, we rewrite ℬ∩i=(xj1(i),yj1(i)),…,(xjKi(i),yjKi(i))B _i=\(x^(i)_j_1,y^(i)_j_1),…,(x^(i)_j_K_i,y^(i)_j_K_i)\, where j1,…,jKi∈[ni]j_1,…,j_K_i∈[n_i] are distinct indices. ⟨,∇Ψ()⟩ π_S,∇ ( π) =⟨∑i=1h1Ki∩i,∇Ψ()⟩=∑i=1h∑k=1Ki1Ki⋅∂Ψ⁡()∂πi,jk=−∑i=1h∑k=1Ki1Ki⋅∂Φ⁡()∂πi,jk = _i=1^h 1K_i1_S _i,∇ ( π) = _i=1^h _k=1^K_i 1K_i· ∂ ( π)∂ _i,j_k=- _i=1^h _k=1^K_i 1K_i· ∂ ( π)∂ _i,j_k (14) =∑i=1h∑k=1Ki⁡[ℒ⁡(∗(−(i,1)))−ℒ⁡(∗(−(i,1)∪(xjk(i),yjk(i))))] = _i=1^h _k=1^K_iE [L ( θ^* (S_-(i,1) ) )-L ( θ^* (S_-(i,1)∪(x^(i)_j_k,y^(i)_j_k) ) ) ] =∑i=1h∑k=1Ki⁡[U⁡(∗(−(i,1)∪(xjk(i),yjk(i))))−U⁡(∗(−(i,1)))] = _i=1^h _k=1^K_iE [U ( θ^* (S_-(i,1)∪(x^(i)_j_k,y^(i)_j_k) ) )-U ( θ^* (S_-(i,1) ) ) ] =∑i=1h∑k=1Ki⁡[U⁡(∗(−(i,k)∪(xjk(i),yjk(i))))−U⁡(∗(−(i,k)))] = _i=1^h _k=1^K_iE [U ( θ^* (S_-(i,k)∪(x^(i)_j_k,y^(i)_j_k) ) )-U ( θ^* (S_-(i,k) ) ) ] where −(i,k)≜⋃(i′,k′)≠(i,k)di′,k′S_-(i,k) _(i ,k )≠(i,k)d_i ,k as the subset obtained by excluding the k-th draw di,kd_i,k from source iD_i and the final inequality comes from the fact that −(i,k)S_-(i,k) and −(i,1)S_-(i,1) are identically distributed. Furthermore, if a set function f is (γ,β)(γ,β)-weakly monotone submodular, we firstly can show that for any two data candidates e1,e2∈e_1,e_2 and the subset ⊆D , we have that γ⁡(f⁡(e1|ℬ∪e2)+f⁡(e2|ℬ))=γ⁡(f⁡(ℬ∪e1,e2)−f⁡(ℬ))≤f⁡(e1|ℬ)+f⁡(e2|ℬ),γ (f(e_1|B∪\e_2\)+f(e_2|B) )=γ (f(B∪\e_1,e_2\)-f(B) )≤ f(e_1|B)+f(e_2|B), such that f⁡(e1|ℬ)≥γf(e1|ℬ∪e2)−(1−γ)f(e2|ℬ).f(e_1|B)≥γ f(e_1|B∪\e_2\)-(1-γ)f(e_2|B). (15) Note that, in the field of set function optimization, we often denote the marginal gain of adding ℬB to A by f⁡(ℬ|)≜f⁡(∪ℬ)−f⁡()f(B|A) f(A )-f(A). Especially when ℬ≜eB \e\ is a singleton, we simply write f⁡(e∣)≜f⁡(e∣)f(e ) f(\e\ ). According to Eq.(15), we subsequently can have that ⟨,∇Ψ()⟩ π_S,∇ ( π) (16) =∑i=1h∑k=1Ki⁡[U⁡(∗(−(i,k)∪(xjk(i),yjk(i))))−U⁡(∗(−(i,k)))] = _i=1^h _k=1^K_iE [U ( θ^* (S_-(i,k)∪(x^(i)_j_k,y^(i)_j_k) ) )-U ( θ^* (S_-(i,k) ) ) ] =∑i=1h∑k=1Ki⁡[U⁡(∗((xjk(i),yjk(i))∣−(i,k)))] = _i=1^h _k=1^K_iE [U ( θ^* ((x^(i)_j_k,y^(i)_j_k) _-(i,k) ) ) ] ≥∑i=1h∑k=1Ki⁡[γ⋅U⁡(∗((xjk(i),yjk(i))∣−(i,k)∪di,k))−(1−γ)⋅U⁡(∗(di,k|−(i,k)))], ≥ _i=1^h _k=1^K_iE [γ· U ( θ^* ((x^(i)_j_k,y^(i)_j_k) _-(i,k)∪ d_i,k ) )-(1-γ)· U ( θ^* (d_i,k|S_-(i,k) ) ) ], where di,kd_i,k is the k-th draw from source iD_i. Then, from the first inequality in Eq.(11), we have that ∑i=1h∑k=1Ki⁡[U⁡(∗((xjk(i),yjk(i))∣−(i,k)∪di,k))] _i=1^h _k=1^K_iE [U ( θ^* ((x^(i)_j_k,y^(i)_j_k) _-(i,k)∪ d_i,k ) ) ] (17) ≥γ⋅[U(∗(∪i=1h∪k=1Ki(x(i)jk,y(i)jk)∣−(i,k)∪di,k))] ≥γ·E [U ( θ^* ( _i=1^h _k=1^K_i(x^(i)_j_k,y^(i)_j_k) _-(i,k)∪ d_i,k ) ) ] =γ⋅[U(∗(ℬ∣∪i=1h∪k=1Kidi,k))] =γ·E [U ( θ^* (B _i=1^h _k=1^K_id_i,k ) ) ] =γ⋅[U(∗(ℬ∪(∪i=1h∪k=1Kidi,k)))−U(∗(∪i=1h∪k=1Kidi,k))] =γ·E [U ( θ^* (B∪( _i=1^h _k=1^K_id_i,k) ) )-U ( θ^* ( _i=1^h _k=1^K_id_i,k ) ) ] ≥γ⋅(U⁡(∗(ℬ))−Ψ⁡()), ≥γ· (U( θ^*(B))- ( π) ), where the final inequality follows from the monotonicity of the utility function U. Furthermore, from the second inequality in Eq.(11), we also have, ∑i=1h∑k=1Ki⁡[U⁡(∗(di,k|−(i,k)))] _i=1^h _k=1^K_iE [U ( θ^* (d_i,k|S_-(i,k) ) ) ] (18) ≤β⋅[U(∗(∪i=1h∪k=1Kidi,k))−U(∗(∅))] ≤β·E [U ( θ^* ( _i=1^h _k=1^K_id_i,k ) )-U( θ^*( )) ] ≤β⋅[U(∗(∪i=1h∪k=1Kidi,k))]=β⋅Ψ(), ≤β·E [U ( θ^* ( _i=1^h _k=1^K_id_i,k ) ) ]=β· ( π), where the second inequality follows from Eq.(11). Merging Eq.(18) and Eq.(17) into Eq.(16), we have that ⟨,∇Ψ()⟩ π_S,∇ ( π) (19) ≥∑i=1h∑k=1Ki⁡[γ⋅U⁡(∗((xjk(i),yjk(i))∣−(i,k)∪di,k))−(1−γ)⋅U⁡(∗(di,k|−(i,k)))] ≥ _i=1^h _k=1^K_iE [γ· U ( θ^* ((x^(i)_j_k,y^(i)_j_k) _-(i,k)∪ d_i,k ) )-(1-γ)· U ( θ^* (d_i,k|S_-(i,k) ) ) ] ≥γ2(U⁡(∗(ℬ))−Ψ⁡())−(1−γ)β⋅Ψ⁡() ≥γ^2 (U( θ^*(B))- ( π) )-(1-γ)β· ( π) =γ2⋅U⁡(∗(ℬ))−(β⁡(1−γ)+γ2)⋅Ψ⁡(). =γ^2· U( θ^*(B))-(β(1-γ)+γ^2)· ( π). Secondly, we investigate ⟨,∇Ψ()⟩ π,∇ ( π) . At first, from the fifth equality in Eq (12) and Eq (13), we have that ∂Ψ⁡()∂πi,j=−∂Φ⁡()∂πi,j ∂ ( π)∂ _i,j=- ∂ ( π)∂ _i,j =−Ki(∑di^,k^∈i^∪∅,(i^,k^)≠(i,1)(∑di,1∈i∪∅ℒ(∗(∪i^,k^di^,k^))∇πi,jp(di,1∣i))∏(i^,k^)≠(i,1)p(di^,k^∣i^)) =-K_i ( _d_ i, k _ i∪\ \,( i, k)≠(i,1) ( _d_i,1 _i∪\ \L ( θ^* ( _ i, kd_ i, k ) ) _ _i,jp(d_i,1 π_i) ) _( i, k)≠(i,1)p(d_ i, k π_ i) ) =−Ki(∑di^,k^∈i^∪∅,(i^,k^)≠(i,1)ℒ(∗((x(i)j,y(i)j)∣∪(i^,k^)≠(k,1)di^,k^))∏(i^,k^)≠(i,1)p(di^,k^∣i^)) =-K_i ( _d_ i, k _ i∪\ \,( i, k)≠(i,1)L ( θ^* ((x^(i)_j,y^(i)_j) _( i, k)≠(k,1)d_ i, k ) ) _( i, k)≠(i,1)p(d_ i, k π_ i) ) =Ki(∑di^,k^∈i^∪∅,(i^,k^)≠(i,1)U(∗((x(i)j,y(i)j)∣∪(i^,k^)≠(k,1)di^,k^))∏(i^,k^)≠(i,1)p(di^,k^∣i^)) =K_i ( _d_ i, k _ i∪\ \,( i, k)≠(i,1)U ( θ^* ((x^(i)_j,y^(i)_j) _( i, k)≠(k,1)d_ i, k ) ) _( i, k)≠(i,1)p(d_ i, k π_ i) ) As a result, we have that ⟨,∇Ψ()⟩ π,∇ ( π) =∑i=1h∑j=1niπi,j⋅∂Ψ⁡()∂πi,j = _i=1^h _j=1^n_i _i,j· ∂ ( π)∂ _i,j =∑i=1h∑j=1niπi,jKi(∑di^,k^∈i^∪∅,(i^,k^)≠(i,1)U(∗((x(i)j,y(i)j)∣∪(i^,k^)≠(k,1)di^,k^))∏(i^,k^)≠(i,1)p(di^,k^∣i^)) = _i=1^h _j=1^n_i _i,jK_i ( _d_ i, k _ i∪\ \,( i, k)≠(i,1)U ( θ^* ((x^(i)_j,y^(i)_j) _( i, k)≠(k,1)d_ i, k ) ) _( i, k)≠(i,1)p(d_ i, k π_ i) ) =∑i=1hKi∑j=1nip(di,1=(x(i)j,y(i)j)∣i)(∑di^,k^∈i^∪∅,(i^,k^)≠(i,1)U(∗((x(i)j,y(i)j)∣∪(i^,k^)≠(k,1)di^,k^))∏(i^,k^)≠(i,1)p(di^,k^∣i^)) = _i=1^hK_i _j=1^n_ip (d_i,1=(x^(i)_j,y^(i)_j) π_i ) ( _d_ i, k _ i∪\ \,( i, k)≠(i,1)U ( θ^* ((x^(i)_j,y^(i)_j) _( i, k)≠(k,1)d_ i, k ) ) _( i, k)≠(i,1)p(d_ i, k π_ i) ) =∑i=1hKi∑j=1ni(∑di^,k^∈i^∪∅,(i^,k^)≠(i,1)U(∗((x(i)j,y(i)j)∣∪(i^,k^)≠(k,1)di^,k^))⋅p(di,1=(x(i)j,y(i)j)∣i)∏(i^,k^)≠(i,1)p(di^,k^∣i^)) = _i=1^hK_i _j=1^n_i ( _d_ i, k _ i∪\ \,( i, k)≠(i,1)U ( θ^* ((x^(i)_j,y^(i)_j) _( i, k)≠(k,1)d_ i, k ) )·p (d_i,1=(x^(i)_j,y^(i)_j) π_i ) _( i, k)≠(i,1)p(d_ i, k π_ i) ) =∑i=1hKi∑j=1ni(∑di^,k^∈i^∪∅,∀i^,k^U(∗(dk,1∣∪(i^,k^)≠(k,1)di^,k^))∏i^=1h∏k^=1Kip(di^,k^∣i^)) = _i=1^hK_i _j=1^n_i ( _d_ i, k _ i∪\ \,∀ i, kU ( θ^* (d_k,1 _( i, k)≠(k,1)d_ i, k ) ) _ i=1^h _ k=1^K_ip(d_ i, k π_ i) ) =∑i=1hKi⋅[U(∗(di,1∣∪(i^,k^)≠(k,1)di^,k^))]=∑i=1h∑k=1Ki⋅[U(∗(di,k∣∪(i^,k^)≠(k,1)di^,k^))] = _i=1^hK_i·E [U ( θ^* (d_i,1 _( i, k)≠(k,1)d_ i, k ) ) ]= _i=1^h _k=1^K_i·E [U ( θ^* (d_i,k _( i, k)≠(k,1)d_ i, k ) ) ] =∑i=1h∑k=1Ki⁡[U⁡(∗(di,k|−(i,k)))]. = _i=1^h _k=1^K_iE [U ( θ^* (d_i,k|S_-(i,k) ) ) ]. Then, from the weak submodularity of U, we have ⟨,∇Ψ()⟩≤∑i=1h∑k=1Ki⁡[U⁡(∗(di,k|−(i,k)))]≤β⋅Ψ⁡(), π,∇ ( π) ≤ _i=1^h _k=1^K_iE [U ( θ^* (d_i,k|S_-(i,k) ) ) ]≤β· ( π), (20) where the final inequality comes from Eq.(11). As a result, by merging Eq. 20 and Eq. 19, we then have that ⟨ℬ−,∇Ψ()⟩ π_B- π,∇ ( π) =⟨ℬ,∇Ψ()⟩−⟨,∇Ψ()⟩ = π_B,∇ ( π) - π,∇ ( π) =≥γ2⋅U(∗(ℬ))−(β(1−γ)+γ2)⋅Ψ()−β⋅Ψ() =≥γ^2· U( θ^*(B))-(β(1-γ)+γ^2)· ( π)-β· ( π) =≥γ2⋅U(∗(ℬ))−(β+β(1−γ)+γ2)⋅Ψ(). =≥γ^2· U( θ^*(B))-(β+β(1-γ)+γ^2)· ( π). ∎ Beside Theorem 4, we next prove the boundedness of ∇Ψ()∇ ( π), that is to say, Lemma 3. if we denote the maximum marginal value of the utility function U(∗())U( θ^*(S)) as MU≜max⊆,e∈∖⁡(U⁡(∗(∪e))−U⁡(∗()))M_U _S ,e (U( θ^*(S∪\e\))-U( θ^*(S)) ), we have that, for any i∈[h]i∈[h], j∈[ni]j∈[n_i], |∂Ψ⁡()∂πi,j|≤MU⋅Kmax, | ∂ ( π)∂ _i,j |≤ M_U· K_max, where Kmax≜maxi∈[h]⁡KiK_max _i∈[h]K_i is the maximum budget over the h communities 1,…,h\D_1,…,D_h\. Similarly, we also can infer that the two-point gradient estimator in Eq.(8) also can be bounded by MU⋅KmaxM_U· K_max. Next, we prove Theorem 3. Proof. At first, let ∗⊆S^* denote the optimal solution for Eq.(6). Then, we denote (t)≜(1(t+1),…,h(t+1)) π^(t) ( π^(t+1)_1,…, π^(t+1)_h) and set the correponding two-point gradient estimator as ∇Φ^() ∇ ( π), namely, ∇Φ^()≜[∂Φ⁡()∂πi,j^]i,j ∇ ( π) [ ∂ ( π)∂ _i,j]_i,j with ∂Φ⁡()∂πi,j^≜Ki(ℒ⁡(∗(−(i,1)∪(xj(i),yj(i))))−ℒ⁡(∗(−(i,1)))) ∂ ( π)∂ _i,j K_i (L ( θ^* (S_-(i,1)∪(x^(i)_j,y^(i)_j) ) )-L ( θ^* (S_-(i,1) ) ) ). Subsequently, according to the projective step in ReSGD, we can show that ‖(t+1)−∗‖22 \| π^(t+1)- π_S^* \|_2^2 ≤‖(t)−η⋅∇Φ^((t))−∗‖22 ≤ \| π^(t)-η· ∇ ( π^(t))- π_S^* \|_2^2 =‖(t)−∗‖22−2η⋅⟨∇Φ^((t)),(t)−∗⟩+η2‖∇Φ^((t))‖22 = \| π^(t)- π_S^* \|_2^2-2η· ∇ ( π^(t)), π^(t)- π_S^* +η^2 \| ∇ ( π^(t)) \|_2^2 =‖(t)−∗‖22−2η⋅⟨∇Φ^((t)),(t)−∗⟩+η2⋅n⋅MU2⋅Kmax2 = \| π^(t)- π_S^* \|_2^2-2η· ∇ ( π^(t)), π^(t)- π_S^* +η^2· n· M_U^2· K_max^2 where the final inequality comes from Lemma 3. As a result, we have that ‖(t+1)−∗‖22≤‖(t)−η⋅∇Φ^((t))−∗‖22 \| π^(t+1)- π_S^* \|_2^2≤ \| π^(t)-η· ∇ ( π^(t))- π_S^* \|_2^2 (21) =‖(t)−∗‖22−2η⋅⁡⟨∇Φ^((t)),(t)−∗⟩+η2⋅n⋅MU2⋅Kmax2 =E \| π^(t)- π_S^* \|_2^2-2η·E ∇ ( π^(t)), π^(t)- π_S^* +η^2· n· M_U^2· K_max^2 =‖(t)−∗‖22−2η⋅⁡[⟨⁡[∇Φ^((t))∣(t)],(t)−∗⟩]+η2⋅n⋅MU2⋅Kmax2 =E \| π^(t)- π_S^* \|_2^2-2η·E [ [ ∇ ( π^(t)) π^(t) ], π^(t)- π_S^* ]+η^2· n· M_U^2· K_max^2 =‖(t)−∗‖22−2η⋅⁡[⟨∇Φ((t)),(t)−∗⟩]+η2⋅n⋅MU2⋅Kmax2 =E \| π^(t)- π_S^* \|_2^2-2η·E [ ∇ ( π^(t)), π^(t)- π_S^* ]+η^2· n· M_U^2· K_max^2 =‖(t)−∗‖22+2η⋅⁡[⟨∇Ψ((t)),(t)−∗⟩]+η2⋅n⋅MU2⋅Kmax2, =E \| π^(t)- π_S^* \|_2^2+2η·E [ ∇ ( π^(t)), π^(t)- π_S^* ]+η^2· n· M_U^2· K_max^2, where the final equality follows from ∇Ψ=−∇Φ∇ =-∇ . Furthermore, according to (γ,β)(γ,β)-weak submodularity, namely, Theorem 4, we also can show that, ‖(t+1)−∗‖22≤‖(t)−∗‖22+2η⋅⁡[⟨∇Ψ((t)),(t)−∗⟩]+η2⋅n⋅MU2⋅Kmax2 \| π^(t+1)- π_S^* \|_2^2 \| π^(t)- π_S^* \|_2^2+2η·E [ ∇ ( π^(t)), π^(t)- π_S^* ]+η^2· n· M_U^2· K_max^2 ≤‖(t)−∗‖22−2η⋅⁡[⟨∇Ψ((t)),∗−(t)⟩]+η2⋅n⋅MU2⋅Kmax2 \| π^(t)- π_S^* \|_2^2-2η·E [ ∇ ( π^(t)), π_S^*- π^(t) ]+η^2· n· M_U^2· K_max^2 ≤‖(t)−∗‖22+2η⋅((β+β⁡(1−γ)+γ2)⋅⁡[Ψ⁡((t))]−γ2⋅U⁡(∗(ℬ)))+η2⋅n⋅MU2⋅Kmax2 \| π^(t)- π_S^* \|_2^2+2η· ((β+β(1-γ)+γ^2)·E[ ( π^(t))]-γ^2· U( θ^*(B)) )+η^2· n· M_U^2· K_max^2 We further can have that ∑t=0T−12η⋅(γ2⋅U⁡(∗(ℬ))−(β+β⁡(1−γ)+γ2)⋅⁡[Ψ⁡((t))]) _t=0^T-12η· (γ^2· U( θ^*(B))-(β+β(1-γ)+γ^2)·E[ ( π^(t))] ) (22) ≤∑t=0T−1(‖(t)−∗‖22−‖(t+1)−∗‖22)+nTη2⋅MU2⋅Kmax2 ≤ _t=0^T-1 (E \| π^(t)- π_S^* \|_2^2-E \| π^(t+1)- π_S^* \|_2^2 )+nTη^2· M_U^2· K_max^2 =‖(0)−∗‖22−‖(T)−∗‖22+nTη2⋅MU2⋅Kmax2. =E \| π^(0)- π_S^* \|_2^2-E \| π^(T)- π_S^* \|_2^2+nTη^2· M_U^2· K_max^2. Due to the randomness of the index τ, we finally can infer that the subset return by our ReSGD satifies that ⁡[U⁡(∗())]=⁡[Ψ⁡((τ))]=∑t=0T−1⁡[Ψ⁡((t))]T [U( θ^*(S))]=E[ ( π^(τ))]= _t=0^T-1E[ ( π^(t))]T ≥(γ2β+β⁡(1−γ)+γ2)⋅U⁡(∗(ℬ))−‖(0)−∗‖22−‖(T)−∗‖222η⋅T⋅(β+β⁡(1−γ)+γ2)−n⋅MU2⋅Kmax22(β+β⁡(1−γ)+γ2)⋅η, ≥ ( γ^2β+β(1-γ)+γ^2 )· U( θ^*(B))- E \| π^(0)- π_S^* \|_2^2-E \| π^(T)- π_S^* \|_2^22η· T·(β+β(1-γ)+γ^2)- n· M_U^2· K_max^22(β+β(1-γ)+γ^2)·η, where the second equality follows from the randomness of the index τ and the final inequality comes from Eq.(22). Since the Lemma B.1 in [102] show that β⁡(1−γ)+γ2≥34β(1-γ)+γ^2≥ 34, when η≜1nKmax2⋅Tη 1nK_max^2· T, we have ⁡[U⁡(∗())]≥γ2β+β⁡(1−γ)+γ2⋅U⁡(∗(∗))−⁡(n⋅Kmax2T)E\! [U( θ^*(S)) ]≥ γ^2β+β(1-γ)+γ^2· U( θ^*(S^*))-O ( n· K_max^2T ). ∎ Appendix E Projection Calculation Algorithm 2 Projection onto Δn _n Input: Vector ∈ℝn z ^n Output: Vector ≜ΠΔn() π _ _n( z) Set +←max⁡, z_+← \ z, 0\; if ⊤+≤1 1 z_+≤ 1 then return + z_+; Find λ∗≥0λ^*≥ 0 such that ∑j=1nmax⁡zj−λ∗,0=1 _j=1^n \z_j-λ^*,0\=1 by bisection. Set πj←max⁡zj−λ∗,0 _j← \z_j-λ^*,0\ for all j∈[n]j∈[n]; return π; In Algorithm 1, each prior vector is updated by projected gradient descent, namely, i(t+1)≜ΠΔni(i(t)−η^i(t)),∀i∈[h], π_i^(t+1) _ _n_i ( π_i^(t)-η g_i^(t) ), ∀ i∈[h], where Δni≜∈[0,1]ni:‖1≤1 _n_i \ p∈[0,1]^n_i:\ \| p\|_1≤ 1\ denotes the sub-probability simplex. If we set ≜i(t)−η^i(t) z π_i^(t)-η g_i^(t), the projection operation is equivalent to solving the following quadratic programming: min∈ℝn⁡12‖−‖22,s.t.≥,⊤≤1. _ π ^n 12\| π- z\|_2^2, .t. π≥ 0, 1 π≤ 1. (23) The main purpose of this section is to provide an efficient method for this singly constrained quadratic programming problem Eq. (23), namely, Algorithm 2. Specifically, Algorithm 2 separates the projection into two cases. It first computes the nonnegative part +≜max⁡, z_+ \ z, 0\. When ‖+‖1≤1\| z_+\|_1≤ 1, this vector already lies in Δn _n and is therefore the Euclidean projection of z onto Δn _n. Specifically, Algorithm 2 proceeds in two cases. It first computes the nonnegative truncation +≜max⁡, z_+ \ z, 0\. If ‖+‖1≤1\| z_+\|_1≤ 1, then + z_+ is already feasible for Eq. (23). In this case, no further normalization is needed, and the projection is simply given by + z_+. If ‖+‖1>1\| z_+\|_1>1, we find the optimal dual variable λ∗≥0λ^*≥ 0 by solving the following one-dimensional equation: ∑j=1nmax⁡zj−λ∗,0=1 _j=1^n \z_j-λ^*,0\=1.Since the left-hand side is monotonically decreasing in λ∗λ^*, this equation can be efficiently solved by bisection [88]. After obtaining λ∗λ^*, the projected vector is computed coordinate-wise as πj≜max⁡zj−λ∗,0,∀j∈[n] _j \z_j-λ^*,0\, ∀ j∈[n]. We now justify Algorithm 2. If +≜max⁡, z_+ \ z, 0\ satisfies ⊤+≤1 1 z_+≤ 1, then +∈Δn z_+∈ _n. Moreover, for each coordinate, the closest nonnegative point to zjz_j is max⁡zj,0 \z_j,0\. Hence, under this case, the nonnegativity projection already satisfies the ℓ1 _1 constraint, and therefore ΠΔn()=+ _ _n( z)= z_+. It remains to consider the case ⊤+>1 1 z_+>1. In this case, the constraint ⊤≤1 1 π≤ 1 must be active at the optimum; otherwise, one could still decrease the objective by moving π closer to + z_+. Therefore, the problem reduces to finding the closest point to z over the set ∈ℝ+n:⊤=1 \ π ^n_+: 1 π=1 \. Consider the following Lagrangian ℒ⁡(,λ,)=12‖−‖22+λ⁡(⊤−1)−⊤,L( π,λ, μ)= 12\| π- z\|_2^2+λ( 1 π-1)- μ π, where λ∈ℝλ is the multiplier for the equality constraint ⊤=1 1 π=1, and ≥ μ≥ 0 is the multiplier for the nonnegativity constraint ≥ π≥ 0. The KKT conditions give πj−zj+λ−μj=0,μjπj=0,πj≥0,μj≥0,∀j∈[n]. _j-z_j+λ- _j=0, _j _j=0, _j≥ 0, _j≥ 0, ∀ j∈[n]. From these conditions, if πj>0 _j>0, then μj=0 _j=0 and thus πj=zj−λ _j=z_j-λ. If πj=0 _j=0, then λ−zj=μj≥0λ-z_j= _j≥ 0 such that zj−λ≤0z_j-λ≤ 0. Therefore, the optimal solution must take the form πj∗=max⁡zj−λ∗,0,∀j∈[n] _j^*= \z_j-λ^*,0\, ∀ j∈[n], for some scalar λ∗λ^*. Since we know that the budget constraint is active in this case, λ∗λ^* is determined by ∑j=1nmax⁡zj−λ∗,0=1 _j=1^n \z_j-λ^*,0\=1. Appendix F Effective Reservoir Rounding Algorithm 3 Reservoir Rounding Scheme Input: Priors ii=1h\ π_i\_i=1^h, budgets Kii=1h\K_i\_i=1^h, datasets ii=1h\D_i\_i=1^h Output: Selected subset ≜⋃i=1h⋃k=1Kidi,kS _i=1^h _k=1^K_id_i,k for i=1,2,…,hi=1,2,…,h do Set Kimin←min⁡‖i‖0,KiK_i ← \\| π_i\|_0,K_i\; Draw KiminK_i distinct samples di,1,…,di,Kimin\d_i,1,…,d_i,K_i \ from iD_i via i π_i without replacement; if Kimin<KiK_i <K_i then Draw Ki−KiminK_i-K_i additional samples di,Kimin+1,…,di,Ki\d_i,K_i +1,…,d_i,K_i\ uniformly at random from i∖di,1,…,di,KiminD_i \d_i,1,…,d_i,K_i \; A straightforward rounding strategy is to independently sample KiK_i data candidates from each source iD_i according to the learned prior i π_i. However, this operation may select the same data candidate multiple times. Consequently, the number of distinct selected data points from source iD_i can be smaller than the prescribed budget KiK_i, i.e., |⋃k=1Kidi,k|<Ki | _k=1^K_i\d_i,k\ |<K_i. When the induced utility function is monotone, such duplicate selections can reduce the practical utility of the final selected subset, since part of the selection budget is effectively wasted. To address this issue, we adopt a simple without-replacement sampling scheme, termed Reservoir Rounding. As detailed in Algorithm 3, Reservoir Rounding first computes Kimin≜min⁡‖i‖0,KiK_i \\| π_i\|_0,K_i\, which is the maximum number of distinct candidates that can be sampled according to the nonzero support of i π_i. It then draws KiminK_i distinct data candidates from iD_i according to i π_i without replacement. If Kimin<KiK_i <K_i, the remaining Ki−KiminK_i-K_i candidates are uniformly sampled from the unselected data points in iD_i. In this way, reservoir rounding scheme always returns exactly KiK_i distinct candidates from each source iD_i. In particular, we show that this reservoir rounding scheme not only ensures that the final selected subset contains exactly KiK_i distinct data candidates from each source iD_i for all i∈[h]i∈[h], but also guarantees that the expected utility of the resulting subset is lossless, that is to say, Theorem 5. Suppose the induced utility function U(∗(ℬ))U( θ^*(B)) is monotone with respect to the selected subset ℬB. For any priors i∈Δni π_i∈ _n_i and ‖i‖1=1\| π_i\|_1=1,∀i∈[h]∀ i∈[h], let S be the subset returned by Algorithm 3. Then S satisfies the following two properties: 1) |∩i|=Ki,∀i∈[h]|S _i|=K_i,∀ i∈[h] and 2) ⁡[U⁡(∗())]≥Ψ⁡(1,…,h)E [U( θ^*(S)) ]≥ ( π_1,…, π_h), where Ψ denotes the probabilistic relaxation of the utility function U(∗(ℬ))U( θ^*(B)). Remark 3. The condition ‖i‖1=1\| π_i\|_1=1 is imposed only for simplicity. When ‖i‖1<1\| π_i\|_1<1, one can naively normalize the prior as ~i≜i‖i‖1 π_i π_i\| π_i\|_1. Since the induced utility function U is monotone, increasing the total sampling mass from ‖i‖1\| π_i\|_1 to 11 can increase the value of the corresponding probabilistic relaxation. Proof. We first verify the exact-budget property. For each source iD_i, Algorithm 3 sets Kimin≜min⁡‖i‖0,KiK_i \\| π_i\|_0,K_i\. It then samples KiminK_i distinct candidates according to i π_i without replacement. If Kimin<KiK_i <K_i, the algorithm additionally samples Ki−KiminK_i-K_i candidates from the remaining unselected data points in iD_i. Therefore, all selected candidates from iD_i are distinct, and |∩i|=Ki,∀i∈[h]|S _i|=K_i,∀ i∈[h]. It remains to prove the utility guarantee. At first, we review the definition Ψ⁡() ( π). Specifically, for each source iD_i, we independently draw KiK_i data candidates ei,1,…,ei,Ki\e_i,1,…,e_i,K_i\ from i∪∅D_i∪\ \. If we let ind≜⋃i=1h⋃k=1Kiei,kS_ind _i=1^h _k=1^K_ie_i,k, we have Ψ⁡()=[U⁡(∗(ind))] ( π)=E_S [U( θ^*(S_ind)) ]. Then, we denote by iA_i the first-stage subset sampled from iD_i without replacement, namely, i≜di,1,…,di,KiminA_i \d_i,1,…,d_i,K_i \. Since the second sampling step only adds extra candidates to ⋃i=1hi _i=1^hA_i, we naturally have that ⁡[U⁡(∗())]≥⁡[U⁡(∗(⋃i=1hi))]E [U( θ^*(S)) ] [U( θ^*( _i=1^hA_i)) ], where S is the output of Algorithm 3. Thus, it suffices to show that ⁡[U⁡(∗(⋃i=1hi))]≥[U⁡(∗(ind))]=Ψ⁡()E [U( θ^*( _i=1^hA_i)) ] _S [U( θ^*(S_ind)) ]= ( π). Specifically, we will first show that, for any subset ~ S independent of the randomness from source iD_i, the following inequality holds i[U⁡(∗(~∪i))]≥⁡[U⁡(∗(~∪(⋃k=1Kiei,k)))],E_A_i [U ( θ^*( S _i) ) ] [U ( θ^*( S∪ ( _k=1^K_ie_i,k )) ) ], (24) where ei,1,…,ei,Ki\e_i,1,…,e_i,K_i\ are drawn independently according to i π_i. If ‖i‖0<Ki\| π_i\|_0<K_i, then the first-stage without-replacement sampling selects all data candidates with nonzero probability under i π_i. Therefore, ⋃k=1Kiei,k⊆i _k=1^K_ie_i,k _i. By monotonicity of U, the above inequality Eq.(24) holds. It remains to consider the case ‖i‖0≥Ki\| π_i\|_0≥ K_i. Let i=pN1(i),…,pNKi(i)A_i=\p^(i)_N_1,…,p^(i)_N_K_i\ where pNk(i)p^(i)_N_k denotes the data pair sampled without replacement from source iD_i. We then prove Eq.(24). It is worth noting that when conditioned on the first Ki−1K_i-1 selected indices N1,…,NKi−1N_1,…,N_K_i-1, the last without-replacement draw selects an unchosen candidate m with probability πi,m1−P,m∉N1,…,NKi−1 _i,m1-P, m∉\N_1,…,N_K_i-1\, where P≜∑r=1Ki−1πi,NrP _r=1^K_i-1 _i,N_r. Hence, [U(∗(~∪viN1,…,viNKi−1∪viNKi))|N1,…,NKi−1] \! [U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\∪\v_i^N_K_i\) ) |N_1,…,N_K_i-1 ] =∑m∉N1,…,NKi−1πi,m1−PU(∗(~∪viN1,…,viNKi−1∪vim)) = _m∉\N_1,…,N_K_i-1\ _i,m1-PU ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\∪\v_i^m\) ) =∑m∉N1,…,NKi−1(πi,m+P⋅πi,m1−P)⋅U⁡(∗(~∪viN1,…,viNKi−1∪vim)) = _m∉\N_1,…,N_K_i-1$ _i,m+ P· _i,m1-P)· U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\∪\v_i^m$ ) ≥∑m∉N1,…,NKi−1πi,m⋅U⁡(∗(~∪viN1,…,viNKi−1∪vim)) ≥ _m∉\N_1,…,N_K_i-1\ _i,m· U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\∪\v_i^m\) ) +P⋅∑m∉N1,…,NKi−1πi,m1−P⋅U(∗(~∪viN1,…,viNKi−1)) + P· _m∉\N_1,…,N_K_i-1\ _i,m1-P· U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\) ) =∑m∉N1,…,NKi−1πi,m⋅U⁡(∗(~∪viN1,…,viNKi−1∪vim))+P⋅U⁡(∗(~∪viN1,…,viNKi−1)) = _m∉\N_1,…,N_K_i-1\ _i,m· U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\∪\v_i^m\) )+P· U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\) ) =∑m∉N1,…,NKi−1πi,m⋅U⁡(∗(~∪viN1,…,viNKi−1∪vim)) = _m∉\N_1,…,N_K_i-1\ _i,m· U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\∪\v_i^m\) ) +∑m∈N1,…,NKi−1πi,m⋅U(∗(~∪viN1,…,viNKi−1∪vim)) + _m∈\N_1,…,N_K_i-1\ _i,m· U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\∪\v_i^m\) ) =ei,Ki[U(∗(~∪viN1,…,viNKi−1∪ei,Ki])), =E_e_i,K_i [U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\∪\e_i,K_i\ ]) ), where the fourth equality follows from the fact that when m∈N1,…,NKi−1m∈\N_1,…,N_K_i-1\, we have U⁡(∗(~∪viN1,…,viNKi−1∪vim))=U⁡(∗(~∪viN1,…,viNKi−1)).U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1∪\v_i^m\\) )=U ( θ^*( S∪\v_i^N_1,…,v_i^N_K_i-1\) ). This shows that the last without-replacement draw can be replaced by an independent draw without increasing the expected function value. Repeating the same argument for the (Ki−1)(K_i-1)-th draw, the (Ki−2)(K_i-2)-th draw, and so on, we obtain the result of Eq.(24). Finally, applying the above source-wise comparison Eq.(24) sequentially over i=h,h−1,…,1i=h,h-1,…,1 and using the tower property of conditional expectation, we obtain ⁡[U⁡(∗(⋃i=1hi))]≥[U⁡(∗(ind))]=Ψ⁡().E [U( θ^*( _i=1^hA_i)) ] _S [U( θ^*(S_ind)) ]= ( π). Combining this with the fact ⁡[U⁡(∗())]≥⁡[U⁡(∗(⋃i=1hi))]E [U( θ^*(S)) ] [U( θ^*( _i=1^hA_i)) ], we have ⁡[U⁡(∗())]≥[U⁡(∗(ind))]=Ψ⁡()E\! [U( θ^*(S)) ] _S [U( θ^*(S_ind)) ]= ( π). ∎ Appendix G Data Summarization on Cifar-100 Dataset (a) (b) (c) (d) Figure 2: Performance comparison between ReSGD and baseline methods on data summarization. (a) Test accuracy across different coreset sizes for ResNet on CIFAR100. (b) and (c) evaluate the robustness of each method against varying rates of pairwise noise and symmetric noise on CIFAR100, respectively. (d) illustrates the performance under different class imbalance factors on CIFAR100. Shaded areas denote the variance across multiple runs. Overall, the results demonstrate the superior performance and stability of ReSGD under various data corruption and distribution shift scenarios. Appendix H Coreset Selection for Existing Continual Learning Methods We instantiate our method by replacing the reservoir sampling component in ER [93] with ReSGD and CS-Rel [104]. Following [104], we evaluate all methods on Split CIFAR-100 and Split Tiny-ImageNet under memory budgets of 200 and 500. Both benchmarks are divided into 10 class-incremental tasks, and all methods use a ResNet-18 backbone [51]. We compare against standard rehearsal-based baselines, including ER, A-GEM [17], and FDR [6]. Results are summarized in Table 2. For Split CIFAR-100, the 100 classes are evenly split into 10 tasks with 10 classes per task. For Split Tiny-ImageNet, the 200 classes are evenly split into 10 tasks with 20 classes per task. In both cases, all samples from the current task are used for training. We adopt the same data augmentation pipeline as Mammoth [12], including random crop, random horizontal flip, and normalization. We also use the same ResNet-18 implementation provided by Mammoth for both datasets. For optimization, we use SGD with batch size 32 on both benchmarks. The learning rate is set to 2×10−22× 10^-2 for Split CIFAR-100 and 3×10−23× 10^-2 for Split Tiny-ImageNet. The cross-entropy loss factor is fixed to 1.0. The number of training epochs is kept consistent with the corresponding baselines. Table 2: Average accuracy (%) on Split CIFAR-100 and Split Tiny-ImageNet. Results are averaged over 10 incremental tasks with memory buffers of 200 and 500 samples. We use ResNet-18 as the backbone for all experiments. Bold indicates the best performance. Methods Split CIFAR-100 Split Tiny ImageNet 200 500 200 500 A-GEM 9.40 ± 0.05 9.42 ± 0.08 8.07 ± 0.08 8.06 ± 0.04 ER 14.18 ± 0.45 21.08 ± 0.16 8.49 ± 0.16 9.99 ± 0.29 FDR 15.32 ± 0.73 22.83 ± 0.73 8.70± 0.19 10.54 ± 0.21 CSReL-ER 15.35 ± 0.73 22.65 ± 0.81 8.66 ± 0.06 10.44 ± 0.17 ReSGD-ER 15.38 ± 0.25 23.21 ± 0.22 8.87 ± 0.05 10.67 ± 0.12 Appendix I Experiment with Adversarial Initialization In this section, we detail the empirical setup designed to rigorously evaluate the error-correction capabilities and robustness of the ReSGD algorithm. The experiment specifically investigates the algorithm’s capacity to filter corrupted labels and recover a high-quality coreset when subjected to an extreme, adversarial initialization state. Dataset Partitioning and Asymmetric Label Corruption. The experiment utilizes a subset of the MNIST dataset, partitioned into two disjoint subsets: Holdout Set holdoutD_holdout consists of 10,00010,000 pristine samples and serves as the reference distribution; Training Pool trainD_train consists of 10,00010,000 separate samples, acting as the active selection pool from which the coreset is sampled. To simulate severe, systematic data corruption, we introduce asymmetric label noise. We uniformly sample a subset ⊂trainC _train containing 1,0001,000 instances. For each instance (xi,yi)∈(x_i,y_i) , the ground-truth label yiy_i is deterministically distorted according to the following mapping: yinoisy=(yiclean+2)(mod10)y_i^noisy=(y_i^clean+2) 10 The indices of the corrupted subset C are recorded for downstream evaluation of the selection purity. Reference Model Pre-training. Prior to the ReSGD sampling loop, a baseline reference model, denoted as θ∗()θ^*(D), is trained to convergence exclusively on the clean dataset holdoutD_holdout. During the subsequent ReSGD iterations, the parameters of θ∗()θ^*(D) are strictly frozen. This model acts as an oracle to compute the reference loss ℓ⁡(xi,yi,θ∗()) (x_i,y_i;θ^*(D)) and does not undergo further parameter updates. Adversarial Probability Initialization and ϵε-Exploration. To stress-test the ReSGD framework, we initialize the sampling probability vector (1)∈ℝ10000 m^(1) ^10000 in a highly adversarial manner. We randomly select a further subset ⊂A , where ||=200|A|=200. To prevent algorithmic stagnation (where clean samples would never be selected or evaluated if assigned a probability of 00), we introduce an ϵε-exploration constant. The initial probability mass is distributed as follows: mi(1)=0.90||,if i∈0.10|train|−||,if i∉m_i^(1)= cases 0.90|A|,&if i \\ 0.10|D_train|-|A|,&if i cases This initialization places 90%90\% of the sampling probability strictly on corrupted samples, forcing the algorithm to recover from a nearly adversarial starting state while maintaining the mathematical viability of gradient computation across the entire pool. Evaluation Metrics To empirically validate the selection efficacy of the proposed method, we monitor three primary metrics: • Distribution Shift: We plot the trajectory of the total probability mass assigned to the corrupted subset C over T iterations: ∑i∈mi(t) _i m_i^(t). We hypothesize a rapid decay from the initial 0.900.90 towards 00. • Coreset Purity: At the final iteration T, we extract the Top-K samples with the highest probability values in (T+1) m^(T+1). The purity is defined as the empirical fraction of clean samples within this Top-K subset. • Downstream Utility: We train a new model from scratch using the final selected Top-K coreset and report its top-1 classification accuracy on the standard MNIST test set, establishing a comparison against a uniform random sampling baseline. Ablation Results. The result is illustrated in Figures 3 and 4. Figure 3: Performance on MNIST under Asymmetric Noise. We track the evolution of the sampling probability mass on corrupted labels (red line) and the coreset purity (blue line) across 50 iterations for various coreset sizes K∈50,…,250K∈\50,…,250\. Despite an adversarial start where 90% of the initial mass is concentrated on noisy samples, ReSGD consistently reallocates probability to clean samples, achieving near-perfect purity across all tested coreset scales. Figure 4: Performance on CIFAR-10 under Asymmetric Noise. Similar to the MNIST experiments, ReSGD demonstrates strong recovery capabilities on the more complex CIFAR-10 dataset. The rapid decay of the noise-associated probability mass (red) and the corresponding climb in the number of clean samples selected (blue) highlight the algorithm’s robustness to systematic label distortion across different target coreset budgets. Appendix J Data Corruption To demonstrate the robustness of selection under noisy conditions, we test data summarization performance under the label noise case on the MNIST, CIFAR-10, and CIFAR-100 datasets. For experiments with label noise and class imbalance, the coreset size is set to 200 for MNIST and CIFAR-10. Here is the detailed protocol: • Label Noise: We experiment with two types of label corruption: Symmetric Noise and Pairwise Noise. For symmetric noise, the true label of a sample is uniformly flipped to any of the other C−1C-1 classes with a probability p. For pairwise noise (pairflip), the true label is flipped to an adjacent class with a probability p. We vary the noise ratio p from 0.2 to 0.8 to test the methods under varying degrees of corruption. • Class Imbalance: We construct a long-tailed (imbalanced) dataset by exponentially decaying the number of samples across classes. Following standard practices, we define the imbalance factor as μ=nmaxnminμ= n_maxn_min, where nmaxn_max and nmin_min represent the number of training samples in the largest (head) and smallest (tail) classes, respectively. Specifically, the number of samples retained for the i-th class is calculated as ni=nmax⋅σin_i=n_max·σ^i, where the decay coefficient is defined as σ=μ−1C−1σ=μ^- 1C-1. In our experiments, we vary the imbalance factor μ from 50 to 200. We analyze a coreset of 200 samples selected within the noise ratio range of 0.0 to 0.4 under pairwise noise, with the results illustrated in Figure 5. Figure 5: Comparison of the number of selected noisy samples across different noise ratios. We evaluate our method against baseline approaches (CSRel [104], Greedy Coreset [8], and PBCS [116]) on three datasets: (a) MNIST using CNN, (b) CIFAR-10 using ResNet, and (c) CIFAR-100 using ResNet. ReSGD consistently selects the fewest noisy samples across all noise levels, demonstrating its superior robustness and ability to distinguish clean data from pairwise noisy samples. Appendix K Experimental Details We note that different baselines in the literature are evaluated under different training and evaluation protocols (e.g., number of training epochs, data repetition strategies, and evaluation schedules). A naive unification of these protocols—for example, adopting the evaluation setting of baseline A for all methods—can lead to suboptimal or misleading performance for other baselines such as B and C, whose effectiveness is tied to their original selection regimes. Therefore, to ensure a fair and faithful comparison, we follow the evaluation protocol used in each original work when reporting results, rather than enforcing a single unified protocol across all methods. To ensure statistical significance, all experimental results are reported as the mean ± standard deviation across 5 independent runs with different random seeds. All experiments were conducted on a server equipped with four NVIDIA H20 GPUs. The total computational cost for reproducing the complete set of experiments is approximately 800 GPU hours. K.1 Data Summarization Experiment Details We conduct CNN-based experiments on MNIST and ResNet-based experiments on CIFAR-10 and CIFAR-100. CNN on MNIST. The CNN architecture comprises two convolutional blocks, each consisting of convolution, dropout, max-pooling, and ReLU activation. The numbers of filters in the two convolutional layers are 32 and 64, respectively, both with a kernel size of (5 × 5). These are followed by two fully connected layers of sizes 128 and 10, with dropout applied. The dropout probability is set to 0.5. We adopt experimental settings consistent with prior work [8, 116, 104]. • Greedy Coreset and Uniform Sampling [8]: The CNN is trained on the data summary using the Adam optimizer with a learning rate of 5e−45e^-4 for 4000 epochs. • PBCS [116]: The projected stochastic gradient descent update for coreset scores is given by t+1←(t−ηℒℬ(∗())∇ln⁡p⁡(∣t)). s^t+1 _C ( s^t- _B ( θ^*( m) ) _ s p ( m s^t ) ). (25) The number of iterations for optimizing the coreset scores s is T=500T=500, where t∈1,2,…,Tt∈\1,2,…,T\, and the learning rate for score updates is η=5e−2η=5e^-2. The inner optimization ∗() θ^*( m) is performed using Adam with a learning rate of 5e−35e^-3 for 300 epochs. The CNN is then trained on the resulting data summary using Adam with a learning rate of 2e−32e^-3 for 4000 epochs. • CSRel [104]: The reducible loss for coreset selection is defined as Gi=ℓ⁡(xi,yi,S)−ℓ⁡(xi,yi,).G_i= (x_i,y_i; θ_S )- (x_i,y_i; θ_D ). (26) The holdout model θ_D is trained using SGD with a learning rate of 1.5e−31.5e^-3 for 125 epochs. The inner model S θ_S is trained using SGD with a learning rate of 2e−22e^-2 for 100 epochs. Finally, the CNN is trained on the selected data summary using SGD with a learning rate of 2e−22e^-2 for 3000 epochs. For ReSGD, the holdout model θ_D is trained using SGD with a learning rate of 1.5e−31.5e^-3 for 200 epochs. The leave-one-out model ∗(−(i,1)) θ^*(S_-(i,1)) is trained using SGD with a learning rate of 5e−45e^-4 for 1000 epochs. Finally, the CNN is trained on the selected data summary using SGD with a learning rate of 2e−22e^-2 for 3000 epochs. The number of iterations T is 20 and η is 1e−71e^-7. ResNet on Cifar-10 and Cifar-100. The experiments of ResNet-18 on Cifar-10 and Cifar-100 follows the setting in [8, 116, 104]. The details of the implementation of baselines are as follows: • Greedy Coreset and Uniform Sampling [8]: The ResNet is trained by Adam optimizer with a learning rate 5e−55e^-5 in 6 epochs. In each epoch, the coreset is repeatedly sampled to 60000 and the total update iterations are 937 (6000/batch size 64). • PBCS [116]: The number of iterations for optimizing the coreset scores s is T=100T=100, where t∈1,2,…,Tt∈\1,2,…,T\, and the learning rate for score updates is η=3e−2η=3e^-2. The inner optimization ∗() θ^*( m) is performed using Adam with a learning rate of 2e−32e^-3 for 80 epochs. The CNN is then trained on the resulting data summary using Adam with a learning rate of 2e−32e^-3 for 4000 epochs. • CSRel [104]: The holdout model θ_D is trained using SGD with a learning rate of 3e−33e^-3 for 100 epochs. The inner model S θ_S is trained using SGD with a learning rate of 1e−21e^-2 for 100 epochs. Finally, the ResNet is trained on the selected data summary using Adam with a learning rate of 5e−45e^-4 for 1800 epochs. Initial coreset size is set to 0 and the selection step is 200. For ReSGD, the holdout model θ_D is trained using SGD with a learning rate of 3e−33e^-3 for 100 epochs. The leave-one-out model ∗(−(i,1)) θ^*(S_-(i,1)) is trained using SGD with a learning rate of 5e−45e^-4 for 100 epochs. Finally, the CNN is trained on the selected data summary using SGD with a learning rate of 5e−45e^-4 for 1800 epochs. The number of iterations T is 20 and η is 1e−81e^-8. K.2 Continual Learning Experiment Details Figure 6: Impact of α on Average Accuracy and Forgetting. Comparison of performance metrics under different loss factors α on the Permuted MNIST dataset. Results show that α∈[0.15,0.20]α∈[0.15,0.20] provides an optimal trade-off for mitigating forgetting while maintaining overall accuracy. We evaluate our data summarization approach for continual learning across four standard benchmarks: Split MNIST [112], Perm MNIST [43], Split CIFAR-10, and Split CIFAR-100. Datasets and Augmentations. Following the evaluation protocol established by [8, 104], the dataset splits are constructed as follows: Split MNIST and Split CIFAR-10 are both divided into 5 tasks with 2 classes per task. Split CIFAR-100 is divided into 10 tasks, each containing 10 classes. Perm MNIST consists of 10 tasks, where each task corresponds to a randomly permuted version of the original MNIST images. For the training data formulation, we randomly sample 1,000 images per task for Split MNIST, Split CIFAR-10, and Perm MNIST. For Split CIFAR-100, we utilize the full dataset for each task. Regarding data augmentation, we apply standard normalization to both Split MNIST and Perm MNIST. For the more complex CIFAR-10 and CIFAR-100 datasets, our augmentation pipeline includes random cropping, random horizontal flipping, and normalization. Network Architectures. We employ distinct backbone architectures tailored to the complexity of each dataset: • Split MNIST: A standard CNN architecture (as detailed in section K.1). • Split CIFAR-10: A standard ResNet-18 architecture (as detailed in section K.1). • Split CIFAR-100: A modified ResNet-18 architecture, wherein a Batch Normalization layer is injected into each convolutional block. • Permuted MNIST: A fully connected network consisting of two hidden layers with 100 units each, ReLU activation functions, and dropout regularization with a probability of 0.20.2 applied to the hidden layers. Optimization and Hyperparameters. Our training procedures follow the methodologies of [8, 116, 104]. The optimization process involves a specific loss factor, denoted as α (refer to the setting in [104]). The hyperparameter configurations are specifically tuned for each dataset as follows: • For Split MNIST, Split CIFAR-10, and Perm MNIST, we utilize the Adam optimizer with a batch size of 256256 and a learning rate of 5×10−45× 10^-4. Models for these three datasets are trained for 400400 epochs. The loss factor α is set to 100.0100.0 for Split MNIST, 20.020.0 for Split CIFAR-10, and 0.150.15 for Perm MNIST. • For Split CIFAR-100, we employ the SGD optimizer with a smaller batch size of 3232 and a higher learning rate of 2×10−22× 10^-2. The model is trained for 100100 epochs, and the loss factor α is set to 4.04.0. This SGD-based setup is consistently applied across all compared methods for the CIFAR-100 benchmark. The choices of the loss factor Following [104], we adopt the forgetting metric to quantify the performance degradation across subsequent tasks for the experiments in Section 5.2. The forgetting for task j after training on task k is computed as: fjk=maxl∈1,…,k−1⁡al,j−ak,j,∀j<k,f_j^k= _l∈1,…,k-1a_l,j-a_k,j, ∀ j<k, (27) where ak,ja_k,j denotes the accuracy on task j upon the completion of task k. We evaluate both average accuracy and forgetting on the Permuted MNIST dataset. As illustrated in Figure 6, we analyze these metrics across varying loss factors α. Empirically, the optimal balance is achieved at α=0.15α=0.15, though α=0.2α=0.2 also yields robust performance in practice. Appendix L Limitations and Broader Impact Limitations. In this work, we introduce a new probabilistic relaxation for the multi-source bilevel data selection problem, which only assigns probability mass to feasible data subsets and thus naturally yields a lossless rounding scheme. In addition, we show that this new probabilistic relaxation admits an effective two-point gradient estimator, which can simultaneously avoid the inverse-Hessian computation required by existing 0−10-1 relaxation and the unbounded score term appearing in prior probabilistic relaxation. Built upon this reducible two-point estimator, we further design a novel stochastic gradient algorithm, termed ReSGD, for the multi-source bilevel data selection problem, and establish a rigorous approximation guarantee for the subset returned by ReSGD under weak submodularity. A limitation of the current analysis is that ReSGD follows the standard random-output scheme in stochastic optimization: it needs to sample τ uniformly from 0,1,…,T−1\0,1,…,T-1\ and returns a subset drawn according to the historical prior (τ) π^(τ). Although theoretically convenient, in practice, this random-output scheme may be inconvenient, since it requires storing all T historical iterates or recovering the sampled historical iterate after training. Fortunately, recent advances [5, 52, 71] in last-iterate convergence of stochastic gradient methods provide a promising way to remove this random-sampling step, so that one can directly return the subset sampled from the final prior (T) π^(T) while hopefully preserving a comparable approximation guarantee. We leave this last-iterate variant of ReSGD as an important direction for future work. Broader Impact. This work mainly focuses on theoretically exploring the multi-source bilevel data selection problem. Therefore, we do not foresee any direct negative societal impact induced by our proposed method.

Instant research discovery