arXiv CS open access

Variable Selection for Feature-Based Newsvendor

2026 · arxiv_cs

arXiv CS · Papers · 2026

deep-learning machine-learning neural-networks

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

Variable Selection for Feature-Based Newsvendor Zhaoliang Yuan School of Artificial Intelligence, The Chinese University of Hong Kong, Shenzhen, 225080014@link.cuhk.edu.cn

Jie Wang School of Artificial Intelligence and School of Data Science, The Chinese University of Hong Kong, Shenzhen, jwang@cuhk.edu.cn

arXiv:2609.01544v1 [stat.ML] 1 Sep 2026

Feature-based newsvendor models use observable covariates to tailor inventory decisions, aiming to balance holding and shortage costs under demand uncertainty. However, high-dimensional feature sets often hinder interpretability and inflate data collection and implementation costs. This paper studies variable selection for the feature-based newsvendor problem under a hard cardinality constraint on the number of selected features. We formulate the resulting ℓ0 -constrained empirical newsvendor problem with ℓ2 -regularization, establish its computational hardness, and develop a mixed-integer second-order cone programming reformulation that strengthens the standard Big-M formulation. To enable scalability beyond exact optimization, we develop a randomized-rounding algorithm with a bi-criteria guarantee and a greedy heuristic. Statistically, we provide theoretical analysis of the resulting sparse policy estimator, including finite-sample estimation error, out-of-sample risk bounds, and support recovery guarantees. Extensive experiments on both synthetic and real data illustrate the computational and statistical trade-offs among various baselines. Our results demonstrate that the proposed variable selection framework achieves competitive out-of-sample operational costs while using substantially fewer covariates. Key Words: Contextual Stochastic Optimization, Variable Selection, Statistical Analysis, Newsvendor

1. Introduction The newsvendor problem is a classic inventory model that seeks to determine the optimal order quantity by balancing the cost of overstocking and understocking in the face of random demand. It is crucial to obtain an accurate estimate of customer demand, which is influenced by a multitude of features, including weather, time, location, and economic conditions. In the big data era, leveraging available data that pairs these influential features (covariates) with corresponding demand observations enables us to make data-driven predictions of future customer demand and more reliable decisions, while ignoring the presence of covariates can lead to inconsistent decisions [3]. A natural paradigm for incorporating features into decision-making is contextual stochastic optimization (CSO) [31]. Broadly speaking, CSO can be divided into three classes: decisionrule optimization [27, 3], estimate-then-optimize [5], and integrated estimation and optimization [14, 21, 30]. In this paper, we focus on the decision-rule optimization approach, aiming to learn an end-to-end mapping from features to decisions. Linear decision rules are 1

attractive because of their simplicity and tractability [3], but they may not capture complex feature-to-decision relationships. To improve out-of-sample performance , nonlinear decision rules based on kernel methods [6, 28] and deep learning [29, 40, 20] have also been studied. Despite this progress, the interpretability of learned decision rules has received comparatively less attention. In many applications, only a small subset of features is truly decisionrelevant. Identifying such features can improve interpretability, reduce data collection and monitoring costs, simplify downstream implementation, and improve robustness when irrelevant features are noisy or unstable. This motivates us to consider variable selection for feature-based newsvendor models. Variable selection has been extensively studied in statistics and machine learning through sparsity-inducing formulations, including Lasso [33], bridge-type penalties [19], the smoothly clipped absolute deviation (SCAD) penalty [37, 16], and the minimax concave penalty (MCP) [39]. These methods can be viewed as tractable surrogates for the ideal but computationally challenging ℓ0 -constrained formulation. Convex penalties are typically easier to optimize, whereas nonconvex penalties may achieve sharper statistical properties but are harder to solve globally. In this paper, we study the following question: How can we efficiently solve variable selection for feature-based newsvendor and what statistical guarantees can be established? Our contributions are summarized as follows. (I) Exact and approximation algorithms. For the variable selection problem, we add an

ℓ2 -norm regularization and reformulate it as a mixed-integer second-order cone program (MISOCP). Compared with the standard Big-M approach for variable selection, our MISOCP formulation yields a tighter continuous relaxation and enables a more efficient algorithm design. Based on the reformulation, we further develop two efficient approximation algorithms and show that the randomized-rounding method admits a bi-criteria guarantee. (II) Statistical guarantees. Under the assumption that demand has sparse and linear dependence on features, we establish the concentration error bound between our estimated policy and the ground truth optimal policy, as well as the out-of-sample performance and support recovery guarantees. These performance bounds break the curse

of dimensionality because they crucially depend on the sparsity level instead of the data dimension. (III) Numerical validation. We compare the exact optimization methods, randomizedrounding, greedy, and Lasso-based approaches using synthetic and retail data, illustrating their computational and statistical trade-offs. 1.1. Related Works Data-Driven Newsvendor. The classical newsvendor model assumes that the demand distribution is known, whereas modern data-driven approaches learn ordering policies directly from historical data. Early work along these lines includes operational-statistics-based approaches [27] and feature-based linear policies for the big-data newsvendor [3]. Subsequent work has developed richer data-driven models based on quantile regression and machine learning [20], reproducing kernel Hilbert space (RKHS)-based policy classes [6, 5], and deep neural networks [29, 40]. Serrano et al. [32] provided a variable selection framework for the feature-based newsvendor problem. They select key features by adding the sparsity budget constraint to the objective and tune penalty parameters using bilevel optimization. Our work contributes to this literature by studying hard-cardinality-constrained policy learning and by providing both computational and statistical guarantees. Chang et al. [9] also studied this framework, but their computational approach can be viewed as a heuristic for solving our ℓ0 -constrained formulation, and their statistical analysis focused on the data privacy viewpoint. Quantile Regression. A fundamental observation in the newsvendor problem is that the optimal order quantity is a conditional quantile of the demand distribution. This connection links contextual newsvendor problems closely to quantile regression, beginning with the seminal work of Koenker and Bassett [23]. Quantile-regression-based methods have been used successfully in feature-based newsvendor models [20]. In high-dimensional settings, sparse quantile regression has been studied under both ℓ1 -regularized and ℓ0 -type formulations [4, 13]. The study in this paper can be interpreted as an ℓ0 -constrained and ℓ2 -regularized empirical quantile regression problem. Variable Selection. Variable selection is a central topic in statistics and machine learning. Classical approaches promote sparsity through penalized formulations [33], or through

greedy and local-search methods [24, 26, 38, 25, 10] that approximately solve ℓ0 -constrained optimization problems. Although ℓ0 -constrained formulations are generally NP-hard and computationally challenging in the worst case, they often perform well in high-dimensional, small-sample, and low-signal-to-noise settings [7, 1, 38]. This motivates the development of scalable algorithms for directly solving the ℓ0 -constrained variable-selection problem for the newsvendor model, as well as theoretical analyses that quantify the optimality gap of approximation algorithms. Compared with classical sparse ridge regression, our problem is further complicated by an empirical loss that is piecewise linear and asymmetric. To address these challenges, we combine perspective reformulations and support-selection algorithms with a statistical analysis tailored to the quantile-regression structure of the newsvendor loss. Notations. For real numbers a, b, we let a ∨ b = max{a, b}. Given a positive integer n, define

[n] = {1, . . . , n}. For a distribution function F (·), we denote F −1 (p) = inf{x ∈ R : F (x) ≥ p}, p ∈ [0, 1]. Given an event E , define 1E (x) = 1 if x ∈ E and otherwise 1E = 0. For a vector P s ∈ Rn , we use s[i] to denote its i-th component, and define its ℓ0 -norm ∥s∥0 = i∈[n] 1s[i]̸=0 . Given an m × n matrix X and a set S ⊆ [n], let XS denote the submatrix of X with columns from the set S . Given a set F , let Conv(F) denote its convex hull, i.e., the smallest convex set containing F .

2. Problem Setup The newsvendor problem seeks to solve the following stochastic optimization problem:

min E[b max{y − z, 0} + h max{z − y, 0}], z≥0

(1)

where the expectation is taken with respect to the random demand y , z denotes the order quantity, and b, h represent the unit back-ordering and holding costs, respectively. A central challenge in solving (1) is that the underlying demand distribution is typically unknown. In the context of big data, auxiliary features can be used for accurate demand estimation. Let x ∈ RD denote a feature vector. For a given x, the optimal order quantity can be expressed as a function f (x) with

n o f (x) = argmin Ey|x b max{y − z, 0} + h max{z − y, 0} . z≥0

Rather than first estimating the conditional demand distribution and then solving the resulting newsvendor problem, one can directly learn an end-to-end policy function f that maps

features to decisions by solving the following contextual stochastic optimization (CSO) problem:

min E(x,y) b max{y − f (x), 0} + h max{f (x) − y, 0} . f ∈F

(2)

The choice of the function class F critically influences both computational tractability and out-of-sample performance. In this paper, we focus on sparse linear decision rules of the form n o F = f (x) = ω ⊤ x + c : ω ∈ RD , c ∈ R, ∥ω∥0 ≤ d . The ℓ0 -constraint ensures that at most d covariates actively affect the decision, thereby facilitating interpretable and explainable policies. Under the linear decision-rule model, the subproblem associated with each fixed support is a convex optimization and therefore enables us to develop mixed-integer convex program reformulations for our variable selection framework. This framework can be extended to nonlinear decision rules through basis expansion. Remark 1 (Linear Demand Model). Consider the linear demand model y = (ω ∗ )⊤ x + ϵ, where ϵ is a random variable independent of x with cumulative distribution function Fϵ . Let τ = b/(b + h). The optimal feature-dependent newsvendor decision is the conditional τ quantile of y given x, namely f ∗ (x) = (ω ∗ )⊤ x + Fε−1 (τ ). Under this linear demand model, the slope in the linear decision rule recovers the true coefficient ω ∗ , while the intercept estimates the τ -quantile of the noise. Remark 2 (Nonnegativity of order quantities). The optimal newsvendor decision has the nonnegative constraint. For the linear decision rule approach, we clip the estimated optimal order quantity to zero if it is negative. When demand is always non-negative, this procedure does not increase the newsvendor loss. In the main analysis, we focus on the unconstrained formulation in (2) for simplicity. The objective in (2) involves the population expectation under the joint distribution of

(x, y). In practice, however, only independent and identically distributed (i.i.d.) samples {(xi , yi )}ni=1 are available. We therefore replace the population expectation with the empirical average and include an ℓ2 -regularization term to improve both computational tractability and out-of-sample performance, which will be discussed in later sections. In summary, it leads to the optimization problem n i 1 1 Xh ⊤ ⊤ b · max{0, yi − ω xi − c} + h · max{0, ω xi + c − yi } + ∥ω∥22 , min F (ω, c) := ω,c: ∥ω∥0 ≤d n i=1 2γ (3)

Problem (3) is nonconvex because of the ℓ0 -constraint, and is therefore generally difficult to solve. In fact, we show that it is NP-hard by a reduction from the exact cover by 3-sets problem. Its proof is provided in Appendix EC.1. Proposition 1 (Computational Hardness). For fixed γ, b, h > 0, Problem (3) is NP-hard whenever the sparsity level d is part of the input. Problem (3) is also related to sparse quantile linear regression. Indeed, it is equivalent to

( min

ω,c: ∥ω∥0 ≤d

) n (b + h) X 1 ρτ yi − (ω ⊤ xi + c) + ∥ω∥22 , n 2γ i=1

b where τ = b+h and ρτ (u) = u τ − 1{u<0} . In other words, the goal of variable selection for the feature-based newsvendor is equivalent to estimating sparse coefficients that minimize the ℓ2 -regularized empirical quantile regression loss. The remainder of the paper is organized as follows. In Section 3, we present a mixedinteger second-order cone reformulation of Problem (3), which can be solved exactly for medium-sized instances. In Section 4, we develop two approximation algorithms with performance guarantees. In Section 5, we establish statistical guarantees for the proposed approach, including consistency, sample complexity, and support recovery.

3. Exact Computational Algorithms In this section, we develop exact algorithms for solving Problem (3) to global optimality. A standard approach is the Big-M reformulation proposed by Dai [13], which converts the sparsity-constrained problem into a mixed-integer convex program. Big-M Method. Recall that the constraint ∥ω∥0 ≤ d requires that at most d entries of ω be nonzero. To model the support of ω , we introduce a binary vector s ∈ {0, 1}D , where s[i] = 1 indicates that ω[i] is allowed to be nonzero. Equivalently, we impose 1{ω[i] ̸= 0} ≤ s[i], i ∈

[D]. Assume there exists a vector M ∈ RD + such that |ω[i]| ≤ M [i], i ∈ [D], for some optimal solution ω . Then the logical constraint can be equivalently enforced by the linear inequality

|ω[i]| ≤ M [i] s[i], i ∈ [D]. This observation leads to the following reformulation of Problem (3).

Proposition 2 (Big-M reformulation [13]). Assume there exists a vector M ∈ RD + such that an optimal solution ω ∗ to Problem (3) satisfies |ω ∗ [i]| ≤ M [i], ∀i ∈ [D]. Then Problem (3) is equivalent to

min ω,c,s

s.t.

n i 1 1 Xh b max{0, yi − ω ⊤ xi − c} + h max{0, ω ⊤ xi + c − yi } + ∥ω∥22 n i=1 2γ X s[i] ≤ d, s ∈ {0, 1}D ,

(Big-M )

i∈[D]

|ω[i]| ≤ M [i] s[i],

∀i ∈ [D].

Big-M method is standard but often computationally expensive for large-scale mixed-integer problems. In practice, such formulations are typically solved by branch-and-bound, which repeatedly solves convex relaxations. A major difficulty is that the Big-M formulation usually yields a weak relaxation: the optimality gap between (Big-M ) and its continuous relaxation, obtained by replacing s ∈ {0, 1}D with s ∈ [0, 1]D , can be substantial. Some heuristics can help derive tight bounds of M and improve numerical performance. For example, let (ω0 , c0 ) be a near-optimal solution to Problem (3), which can be obtained using some tractable algorithms outlined in Section 4. Any optimal solution to Problem (3) satisfies that

∥ω ∗ ∥22 ≤ F (ω ∗ , c∗ ) ≤ F (ω0 , c0 ). 2γ p Based on this relation, the vector M with M [i] = 2γF (ω0 , c0 ), i ∈ [D] is a valid choice because |ω ∗ [i]| ≤ ∥ω ∗ ∥2 ≤ M [i]. However, solving the Big-M formulation remains challenging for large-scale datasets. Motivated by this limitation, we next introduce an alternative reformulation that is more computationally tractable. We show that it has a tighter convex relaxation than the Big-M approach and is therefore better suited for exact computation. Perspective Reformulation. We introduce a stronger reformulation than the Big-M approach, based on the perspective function. To this end, we first rewrite the quadratic regularization term by introducing auxiliary variables µ ∈ RD +:

1 ∥ω∥22 = min 2γ µ∈RD +

∥µ∥1 2 : ω[i] ≤ µ[i], ∀i ∈ [D] . 2γ

Therefore, Problem (3) is equivalent to the following extended Big-M formulation:

min

ω,c,s,µ

s.t.

n i ∥µ∥ 1 Xh 1 b · max{0, yi − ω ⊤ xi − c} + h · max{0, ω ⊤ xi + c − yi } + , n i=1 2γ X s[j] ≤ d, s ∈ {0, 1}D ,

(4a) (4b)

j∈[D]

(ω[j], s[j], µ[j]) ∈ Sj ,

(4c)

∀j ∈ [D],

where for each coordinate j ∈ [D], define the set

n o 2 Sj := (ω[j], s[j], µ[j]) ∈ R × {0, 1} × R+ : ω[j] ≤ µ[j], |ω[j]| ≤ M [j] s[j] . Günlük and Linderoth [18] have characterized the convex hull of Sj as

n o Conv(Sj ) = (ω[j], s[j], µ[j]) ∈ R × [0, 1] × R+ : ω[j]2 ≤ µ[j]s[j], |ω[j]| ≤ M [j] s[j] . By convexifying the constraint (4c), we can obtain the perspective reformulation of (3). Proposition 3 (Perspective Reformulation). Under the same condition as in Proposition 2, Problem (3) is equivalent to

min

ω,c,s,µ

s.t.

n i 1 X 1 Xh b · max{0, yi − ω ⊤ xi − c} + h · max{0, ω ⊤ xi + c − yi } + µ[j], n i=1 2γ j∈[D] X s[j] ≤ d, s ∈ {0, 1}D , j∈[D]

(ω[j], s[j], µ[j]) ∈ Conv(Sj ),

∀j ∈ [D]. (Perspective)

Since the constraint ω[j]2 ≤ µ[j] s[j] is representable by a rotated second-order cone when

s[j] ≥ 0 and µ[j] ≥ 0, Problem (Perspective) is a mixed-integer second-order cone program (MISOCP). The following theorem compares the Big-M and perspective formulations. Its proof follows techniques in [38] and is provided in Appendix EC.2. Theorem 1 (Strength of Perspective Reformulation).

(I) The mixed-integer formula-

tions (4) and (Perspective) are equivalent, and hence have the same optimal value; Big-M Pers (II) Let Vrel and Vrel denote the optimal values of the continuous relaxations of (Big-M ) Big-M Pers and (Perspective), respectively. It always holds that Vrel ≤ Vrel .

Remark 3 (Perspective Reformulation is Strictly Stronger). We provide an example to show that the perspective reformulation is strictly stronger. Take D = 2, d =

1, b = h = 1, n = 4, M = 2, γ = 1. Assume that one has access to four data samples: {(e1 , 1), (−e1 , −1), (e2 , 1), (−e2 , −1)}, where e1 , e2 denote the unit vectors in R2 . In this case, Big-M Pers it is easy to verify Vrel = 0.75 and Vrel = 0.875. Hence, this problem instance ensures that Big-M Pers Vrel < Vrel .

In summary, perspective reformulation is preferred to the Big-M method in two ways. First, Theorem 1 together with Remark 3 explains that the perspective reformulation is more computationally efficient. Since its continuous relaxation is tighter than that of the Big-M formulation, a branch-and-bound solver typically explores fewer nodes and converges faster in practice. Second, it is difficult to derive valid and tight Big-M constants for solving both formulations, and using overly conservative Big-M constants makes the formulations numerically ill-conditioned to solve. Even if we drop the constraints |ω[j]| ≤ M [j]s[j], j ∈ [D] in Problem (Perspective), it is still a valid reformulation of Problem (3) and computationally tractable, leading to the following Big-M-free reformulation: n i 1 X 1 Xh b max{0, yi − ω ⊤ xi − c} + h max{0, ω ⊤ xi + c − yi } + µ[j] ω,c,s,µ n i=1 2γ j∈[D] X s.t. s[j] ≤ d, s ∈ {0, 1}D ,

min

j∈[D]

ω[i]2 ≤ µ[i]s[i], i ∈ [D],

µ ∈ RD +.

(Big-M -free-Perspective)

4. Scalable Approximation Algorithms In this section, we study two scalable approximation algorithms to obtain near-optimal solutions to Problem (3) and build their performance guarantees. 4.1. Continuous Relaxation with Randomized Rounding Our first approach is motivated by the MISOCP (Big-M -free-Perspective). Its continuous relaxation replacing the binary constraint s ∈ {0, 1}D with s ∈ [0, 1]D becomes a convex program, and can be efficiently solved using off-the-shelf solvers. Then, we estimate the support of sparse regression coefficients using a randomized rounding scheme. See Algorithm 1 for detailed implementation.

Algorithm 1 Continuous Relaxation with Randomized Rounding for Problem (3) Input: Data {(xi , yi )}ni=1 , sparsity budget d, regularization parameter γ 1: Solve the continuous relaxation of (Big-M -free-Perspective) and obtain an optimal solu-

tion (b ω, b c, sb, µ b) 2: Initialize set S = ∅ 3: for i = 1, . . . , D do 4:

Sample an uniform number U on [0, 1].

Update S = S ∪ {i} if U ≤ sb[i].

6: end for 7: Obtain (ω̃, c̃) as the optimal solution to arg min ω,c

8: return

n o F (ω, c) : ω[i] = 0 if i ∈ /S .

(ω̃, c̃) as the regression coefficient estimator

We provide performance guarantees for this computationally efficient algorithm below. Let X = (x1 , . . . , xn )⊤ ∈ Rn×D be the data matrix. For k ∈ [D], define

θk := max σmax (XS X⊤ S ), |S|=k

where σmax (·) denotes the maximum singular value. Due to the randomized rounding procedure, the support size of the estimated regression coefficient may not satisfy the sparsity budget. The following theorem provides the bi-criteria approximation guarantee of Algorithm 1. It proves that with high probability, the selected support is only moderately larger than d, while the objective value is within an additive tolerance of the optimal value. Theorem 2 (Approximation Guarantee of Algorithm 1). Given δ ∈ (0, 1) and ϵ > 0. Let

(ω̃, c̃) be the output from Algorithm 1. Suppose the regularization coefficient p log(2n/δ)(h ∨ b)2 θ1 (h ∨ b)2 2θ1 θd log(2n/δ) 1 ≥ + . 2γ 6nϵ 4nϵ Then, with probability at least 1 − δ , it holds that ! r 2 log(2/δ) log(2/δ) ∥ω̃∥0 ≤ 1 + + d, d d

F (ω̃, c̃) ≤ optval(3) + ϵ.

The proof of Theorem 2 is provided in Appendix EC.3. The proof step to bound the sparsity level of ω̃ follows from the usage of concentration inequality discussed in [38, Theorem 5]. The reference therein implicitly assumed that log(2/δ) ≤ d/3. Our concentration result refines

the concentration result and does not impose this assumption. The support-size bound in Theorem 2 depends on the target sparsity d, but not on the data dimension D. The objec1 tive performance guarantee requires that the regularization parameter 2γ is not too small: its

lower bound is related to the dimension D through the maximum singular values θ1 and θd from the design matrix. This assumption is not quite restrictive since adding ℓ2 -norm regularization ensures the estimated solution is more robust than unregularized models [8]. 4.2. Greedy Algorithm Greedy algorithm is a common heuristic for many variable selection problems. To develop a greedy algorithm for solving (3), we rewrite it as a combinatorial formulation. Indeed, by introducing a subset S ⊆ [D], Problem (3) is equivalent to min min F (ω, c) : ω[i] = 0 if i ∈ /S . S⊆[D],|S|≤d

(5)

ω∈RD ,c∈R

For a fixed support S , the inner subproblem of (5) is a convex program. Using strong duality theory on the inner minimization problem, this subproblem has an equivalent dual representation that depends on S only through submatrix XS . Theorem 3 (Minimax Reformulation of (5)). Problem (5) has the same optimal value as the following minimax optimization problem: min G(S) := maxPn S⊆[D],|S|≤d

z[i]∈[−b,h],i∈[n],

−z

i=1 z[i]=0

⊤

γ 2n

X X⊤ 2 S S

1 ⊤ z− z y . n

(6)

Based on Theorem 3, Greedy algorithm proceeds by iteratively finding a candidate index such that the objective G is minimized for current subset including the new index. See the detailed procedure in Algorithm 2. Algorithm 2 Greedy Algorithm for Problem (6) Input: Data {(xi , yi )}ni=1 , sparsity budget d, regularization parameter γ 1: Initialize S0 = ∅ 2: for k = 0, 1, . . . , d − 1 do 3:

Choose jk ∈ arg minj∈[D]\Sk G(Sk ∪ {j})

Update Sk+1 = Sk ∪ {jk }

5: end for 6: Obtain (ω̃ gr , c̃gr ) as the optimal solution to arg min ω,c

7: return

(ω̃ gr , c̃gr )

n o F (ω, c) : ω[i] = 0 if i ∈ / Sd .

By Theorem 3, it can be shown that the set function G(S) is non-increasing in S , and the objective of iterates in Algorithm 2 is monotonically decreasing. It is an open question on the approximation guarantees of the final iterate from the greedy algorithm, as the set function

G does not satisfy the submodular assumption commonly used in literature. Inspired by [38], the greedy approach can also be integrated with Algorithm 1. Specifically, given threshold

δ > 0 and ω b obtained from continuous relaxation, define C = {i : |b ω [i]| > δ} as the effective support obtained from continuous relaxation. We can apply the greedy algorithm based on C instead of [D] and can reduce the number of candidate variables. This algorithm works well especially when ω b is sparse.

5. Statistical Performance Guarantees In this section, we investigate statistical performance guarantees for our proposed variable selection framework. Assume the observed data samples {(xi , yi )}ni=1 are i.i.d. copies of

(X, Y ), where Y = (ω ∗ )⊤ X + c∗ + ε.

(7)

Here, ω ∗ ∈ RD is a fixed d-sparse vector with S ∗ = supp(ω ∗ ), and |S ∗ | = d, and c∗ ∈ R is a fixed intercept parameter. The distribution of the noise ε is independent of X . We analyze the estimator from the optimization problem ( n ) i ∥ω∥2 1 Xh 2 (b ωn , b cn ) = arg min b max{0, yi − ω ⊤ xi − c} + h max{0, ω ⊤ xi + c − yi } + , n 2γ n ω,c: ∥ω∥0 ≤d i=1 where the regularization parameter γn is allowed to vary with the sample size n. We impose the following technical assumptions throughout our analysis. Assumption 1 (Structure of Data).

(I) The distribution of ε satisfies Fϵ (0) := P(ϵ ≤ 0) =

τ , where τ = b/(b + h). (II) The random variable ε satisfies E|ε| < ∞. (III) The density of the noise exists and is denoted as fε (·). There exists r0 such that for all

|u| ≤ r0 , it holds that fε (u) ≥ f > 0. e = (X ⊤ , 1)⊤ is bounded almost surely: ∥X∥ e 2 ≤ K. (IV) The augmented covariate X (V) There exists κ > 0 such that for each u ∈ RD , t ∈ R, ∥u∥0 ≤ 2d, it holds that

E (u⊤ X + t)2 ≥ κ ∥u∥22 + t2 .

It is noteworthy that Assumption 1(I) is not restrictive. Indeed, (7) can be reformulated as

Y = (ω ∗ )⊤ X + c̃∗ + ε̃∗ for constant c̃∗ = c∗ + Fε−1 (τ ) and random noise ε̃∗ = ε − Fε−1 (τ ). Such a scalar transformation always ensures Assumption 1(I) holds. Assumption 1(II) is a sufficient condition that ensures the optimal newsvendor risk E(X,Y ) b max(Y − (ω ∗ )⊤ X − c∗ , 0) + h max((ω ∗ )⊤ X + c∗ − Y, 0) = (b + h)Eε [ρτ (ε)] is finite. Assumption 1(III) imposes local regularity of the noise density around the target quantile. Assumption 1(IV) imposes the boundedness of the feature vector. Assumption 1(V) is a restricted eigenvalue condition, which ensures identifiability over sparse directions. These assumptions are commonly used in sparse statistical analysis [35]. We first establish a finite-sample estimation-error bound regarding the estimator (b ωn , b cn ) in Theorem 4. It reveals that under proper scaling of the ridge parameter, the estimator converges to the ground truth with a rate crucially dependent on the sparsity level instead of the data dimension. The estimation error vanishes quickly as the sample size n increases. Theorem 4 (Estimation Error Bound). Let δ ∈ (0, 1) and suppose Assumption 1 holds. Set q the ridge parameter 2γ1n ≲ d log(eD/d) , then with probability at least 1 − δ , it holds that n   ! r ω bn − ω ∗ d log(eD/d) + log(1/δ)  =O  , n b cn − c∗ 2

where O(·) hides the constant depending only on b, h, r0 , f , κ, K, ∥ω ∗ ∥2 . The proof of Theorem 4 is provided in Appendix EC.4. Building on Theorem 4, we further show the out-of-sample performance and support recovery guarantees of our obtained estimator. Given regression coefficients (ω, c), define its population newsvendor risk R(ω, c) = E(X,Y ) b max(Y − ω ⊤ X − c, 0) + h max(ω ⊤ X + c − Y, 0) .

(8)

Our out-of-sample performance guarantee in Corollary 1 quantifies the gap between the population newsvendor risk of (ω̂n , ĉn ) and the optimal risk, which vanishes in the order of

d log(eD/d)/n. Corollary 1 (Out-of-Sample Performance Guarantee). Under the same setup as in Theorem 4, with probability at least 1 − δ , it holds that

0 ≤ R(ω̂n , ĉn ) − min

∥ω∥0 ≤d,c

R(ω, c) ≤ O

d log(eD/d) + log(1/δ) , n

where O(·) hides the constant depending only on b, h, r0 , f , κ, K, ∥ω ∗ ∥2 .

The support recovery guarantee is provided in Corollary 2, whose proof is provided in Appendix EC.6 and uses the fact that ω bn converges to the true sparse vector ω ∗ , and so is its support. Corollary 2 (Support Recovery Guarantee). Under the same setup as in Theorem 4, and assume the sample size n is sufficiently large so that with probability at least 1 − δ , it holds that   ∗ ω b −ω  n  < min |ω ∗ [j]|. ∗ j∈S ∗ b cn − c 2

Then, with probability at least 1 − δ ,

ω bn [i] ̸= 0,

ω bn [j] = 0,

∀i ∈ S ∗ , j ∈ / S ∗.

6. Numerical Study In this section, we evaluate the performance of the proposed methods. We generate synthetic data from the linear model Y = (ω ∗ )⊤ X + ε, where ε ∼ N (5, σ 2 ). The feature vector X follows a multivariate Gaussian distribution with mean zero and covariance matrix Σ, where

Σ[i, j] = 0.5|i−j| , i, j ∈ [D]. Although our statistical analysis assumed that the feature vector has bounded support, we use Gaussian-distributed feature vector with unbounded support for numerical evaluation. We further perform scalar transformation to ensure Assumption 1(I) ,i ∈ holds. The first d entries of ω ∗ are nonzero and are specified as ω ∗ [i] = (−1)i+1 1 − i−1 d

[d]. We set the unit shortage and holding costs to b = 2 and h = 1, respectively. For baseline comparison, our experiment includes the following methods: (I) the exact algorithms based on the Big-M and perspective formulations in Section 3; (II) the continuousrelaxation algorithm described in Algorithm 1; (III) the greedy algorithm described in Algorithm 2; and (IV) an ℓ1 -based baseline obtained by replacing the ℓ0 -constraint in Problem (3) with an ℓ1 -regularized formulation. For each algorithm, we impose a time limit of 1800 seconds. The exact algorithm baseline is built on the perspective formulation except in Section 6.1. The mixed-integer optimization algorithms terminate once the relative optimality gap falls below 1e-4. Unless otherwise specified, the ridge and Lasso tuning parameters are selected by cross-validation. We measure the performance support recovery using false discovery proportion (FDP) and non-discovery proportion (NDP), defined in [2]:

|S \ S ∗ | FDP(S) = , |S|

|S ∗ \ S| NDP(S) = , |S ∗ |

(9)

where S ∗ denotes the true support and S denotes the support estimated by a given algorithm. Smaller values of FDP and NDP indicate better support recovery performance. All computations are performed in 2.60GHz Intel(R) Xeon(R) Platinum 8358 CPU and 32GB main memory with Gurobi 13.0.1 solver. 6.1. Comparison of exact algorithms We first compare the two exact algorithms in Section 3 under different configurations of

(D, n, d). Throughout this subsection, we fix γ = 50, M = 100, and σ 2 = 1. Performance is evaluated in terms of the objective value, computational time, and FDP. Table 1 reports the computational results for solving (Big-M ) and (Perspective) via branch-and-bound. Overall, the perspective formulation outperforms the Big-M formulation on most instances, and the advantage becomes more pronounced as the ambient dimension D increases. This observation is consistent with Theorem 1, which shows that (Perspective) provides a stronger continuous relaxation. When D ∈ {1000, 2000}, both approaches often reach the prescribed time limit, indicating that solving the sparse newsvendor problem to global optimality can be computationally demanding in large-scale settings. We also observe that the FDP generally decreases as

D increases in these experiments. A possible explanation is that the regularization parameter γ is not optimally tuned such that the FDP values for different instances are not comparable. Table 1

Computational results of exact algorithms. Here we fix cardinality d = 50. NA indicates the algorithm fails to find any feasible solution within time limit.

Big-M

Perspective

Objective Gap (%)

Time (s)

FDP (%)

Objective Gap (%)

Time (s)

FDP (%)

100

500 1000 5000

10.71 10.52 11.06

3.77 1.02 0.38

1800.00 1800.00 1800.00

52.00 48.00 26.00

10.65 10.45 10.91

9.69e-3 9.94e-3 0.14

11.33 50.79 1800.00

42.00 40.00 24.00

500

500 1000 5000

12.38 10.48 10.77

100 36.71 6.85

1800.00 1800.00 1800.00

12.00 11.00 7.50

10.01 10.11 10.71

0.56 1.37 1.39

1800.00 1800.00 1800.00

9.70 9.70 5.20

500 1000 1000 5000

9.55 11.54 10.81

100.00 100.00 15.52

1800.00 1800.00 1800.00

5.15 7.63 4.26

9.17 10.34 10.74

1.44 1.44 3.17

1800.00 1800.00 1800.00

4.85 4.15 3.68

500 2000 1000 5000

12.01 12.98 NA

100 100 NA

1800.00 1800.00 1800.00

6.41 5.51 NA

9.90 9.56 11.37

1.99 4.34 41.04

1800.00 1800.00 1800.00

2.46 2.41 2.46

6.2. Results on Support Recovery Next, we compare the variable-selection performance of all methods. The noise variance σ 2 ∗ ⊤

) X) is chosen so that Var((ω = 10, that is, the signal-to-noise ratio (SNR) is 10. Var(ε)

D=100

D=300

D=500

1.00

0.75

Methods 0.75

Methods Exact Algorithm

FDP

0.50

Continuous Relaxation 0.50 Greedy Algorithm

0.25

1000

2000

3000

5000

1000

2000

3000

5000

500

Sample Size n

2000

3000

5000

Sample Size n

n=5000

n=1000 1.00

1.00

0.75

Methods 0.75

Exact Algorithm

Continuous Relaxation 0.50

Continuous Relaxation

0.25

Greedy Algorithm

0.25 Lasso

Lasso

0.00

FDP

0.50

Continuous Relaxation 0.50

0.00

0.00 100

200

300

Data Dimension D

500

Methods

Exact Algorithm

FDP

1000

FDP with different choices of sample size n and each subplot represents choices of dimension D.

n=500

Figure 2

Greedy Algorithm Lasso

1.00

Continuous Relaxation

0.00

500

Sample Size n

Figure 1

Continuous Relaxation 0.50

0.25 Lasso

0.00 500

Exact Algorithm

Greedy Algorithm

0.25 Lasso

0.00

Methods

0.75

Exact Algorithm

100

200

300

Data Dimension D

500

100

200

300

500

Data Dimension D

FDP with different choices of dimension D and each subplot represents choices of sample size n.

We first consider the case in which the sparsity budget matches the true support size and is fixed at 10. In this setting, FDP and NDP coincide, so we report only FDP. Figures 1 and 2 show the FDP values under different combinations of (n, D). The exact algorithms achieve the best support recovery performance, especially in the high-dimensional and small-sample regime. Among the scalable baselines, the greedy algorithm generally outperforms both the continuous-relaxation and Lasso baselines.

Plot of FDP (𝑛 = 500)

Figure 3

Plot of NDP (𝑛 = 500)

Plot of FDP (𝑛 = 2000)

Plot of NDP (𝑛 = 2000)

Comparison in terms of FDP and NDP for n ∈ {500, 2000} and D = 100.

We next fix the true support size at 10 and vary the sparsity budget to compare the methods in terms of both FDP and NDP. Figure 3 shows the results for (D, n) ∈

{(100, 500), (100, 2000)}. The black dashed line represents the ground-truth optimum of these metrics. As shown in the figure, exact algorithm remains close to the oracle benchmark across a wide range of sparsity budgets. 6.3. Results on Parameter Estimation We examine the performance of parameter estimation for various methods on problem instances where D = 100, n ∈ {k · 100, k ∈ [10]}, and SNR ∈ {0.1, 1, 10}. The ℓ1 or ℓ2 regularization parameters are tuned using cross-validation. Figure 4 summarizes the numerical results, with the sample size n on the x-axis and the logarithm of relative error on the y axis. The relative error is defined as |ĉ − c∗ |/|c∗ | or ∥ω̂ − ω ∗ ∥2 /∥ω ∗ ∥2 , where (ω̂, ĉ) denotes the estimated parameters and (ω ∗ , c∗ ) the ground truth. All box plots are generated using 5 independent trials. The results indicate that, for all methods, estimation accuracy improves as the sample size n increases. In particular, the exact and greedy algorithms consistently outperform the continuous relaxation and Lasso-based approaches in estimating the weight vector ω ∗ , especially in the small-sample and low-SNR regimes. 6.4. Results on Out-of-Sample Performance We examine the out-of-sample performance of various methods using the population newsvendor risk defined in (8), which is accurately estimated via sample average approxi-

Exact Algorithm

Continuous Relaxation

Greedy Algorithm

SNR = 0:1, D = 100, c ¤ = 6:684111

SNR = 0:1, D = 100

-2

1.5

log

j^c ¡ c ¤ j jc ¤ j

jj! ¡ ! ¤ jj 2 jj! ¤ jj 2

-4

log

Lasso

-6

1.0 0.5 0.0

-8

-0.5 100

200

300

400

500 600 700 Sample size n

800

900

1000

100

200

300

SNR = 1:0, D = 100, c ¤ = 5:532563

400

500 600 700 Sample size n

800

900

1000

800

900

1000

800

900

1000

SNR = 1:0, D = 100 0.0

-3 -0.5

j^c ¡ c ¤ j jc ¤ j

jj! ¡ ! ¤ jj 2 jj! ¤ jj 2

-4

-1.0

log

-5 -6

-1.5

-7 -8

-2.0 100

200

300

400

500 600 700 Sample size n

800

900

1000

100

SNR = 10:0, D = 100, c ¤ = 5:168411

500 600 700 Sample size n

SNR = 10:0, D = 100

jj! ¡ ! ¤ jj 2 jj! ¤ jj 2

-1.5

-6

-2.0

-7

log

j^c ¡ c ¤ j jc ¤ j

log

400

-1.0

-5

-8

-2.5 -3.0

-9

Figure 4

300

-0.5

-4

-10

200

-3.5 100

200

300

400

500 600 700 Sample size n

800

900

1000

100

200

300

400

500 600 700 Sample size n

Comparison in terms of logarithm of relative errors for different (n, D, SNR).

mation with 105 testing samples. The experimental setup follows that of the previous subsection, except that we also examine different choices of the sparsity level d. Recall the ground truth weight parameter ω ∗ is fixed with support size |S ∗ | = 10. For the Lasso method, the

ℓ1 regularization parameter is fine-tuned to achieve the desired sparsity level in the estimated parameter. Figure 5 summarizes the numerical results, with all error bars generated from 5 independent trials. The left-hand-side subplots display the sample size n on the xaxis. These results show that all methods improve with increasing n, while the exact and greedy algorithms consistently yield lower out-of-sample operational costs compared to the continuous relaxation and Lasso-based approaches. The right-hand-side subplots present the

sparsity level d on the x-axis. When the true sparsity level |S ∗ | is underestimated (d < |S ∗ |), the exact and greedy algorithms dominate the others in terms of performance. When it is overestimated, the out-of-sample performance tends to stabilize. Together with the results in Sections 6.2 and 6.3, our numerical findings demonstrate that the variable selection framework delivers satisfactory parameter estimation and out-of-sample performance while using substantially fewer covariates and yielding interpretable decisions.

D = 100, SNR = 0:1 12

4 100

200

300

400

500 600 Sample size n

700

800

900

1000

1.6

1.4 100

200

300

400

500 600 Sample size n

700

800

900

D = 100, SNR = 10:0 0.9 0.8 0.7 0.6 0.5 100

Figure 5

200

300

400

500 600 Sample size n

700

800

900

3.875 3.850 3.825 3.800 3.775 3.750

1000

1.8

Lasso

3.900

1000

2.0

Greedy Algorithm

ntrain = 3000, D = 100, SNR = 0:1

Average out-of-sample operational cost

D = 100, SNR = 1:0

Average out-of-sample operational cost

Continuous Relaxation

Average out-of-sample operational cost

Exact Algorithm

9 10 11 12 13 14 15 16 17 18 19 20 Sparsity level d

ntrain = 3000, D = 100, SNR = 1:0 1.9 1.8 1.7 1.6 1.5 1.4 1

9 10 11 12 13 14 15 16 17 18 19 20 Sparsity level d

ntrain = 3000, D = 100, SNR = 10:0 1.4 1.2 1.0 0.8 0.6 0.4 1

9 10 11 12 13 14 15 16 17 18 19 20 Sparsity level d

Comparison in terms of out-of-sample operational cost for different (n, d, SNR).

7. Real-world case study In this section, we evaluate the performance of various approaches using retail sales data from Kaggle[17]. The real-world dataset is provided by Corporación Favorita, an Ecuadorianbased grocery retailer that sells over 200,000 products. During preprocessing, we focus on hardware products and define the response variable Y as the sum of weekly unit sales across all four hardware product types. The covariate vector X is constructed using standard feature transformations, including one-hot encoding and normalization. The data are restricted to the time period from 2012-12-31 to 2013-12-30. The raw features in the dataset include: store number, city index, week index, month index, weekly averaged oil price, weekly promotion intensity, weekly transaction volume, holiday indicator, weekly number of hardware SKUs sold, and weekly number of active sales days. In addition, we incorporate the following two-way interactions as suggested in Chang et al. [9]: (i) holiday indicator with weekly promotion intensity; (ii) store number with weekly averaged oil price; (iii) month index with city index; (iv) month index with store number; and (v) holiday indicator with store number. After preprocessing, the final dataset consists of n = 2239 samples and D = 1157 dimensions. We split the data into 70% training and 30% testing sets across 20 independent trials. The training data are used to estimate the weight vector and intercept, while the testing data are used for evaluating the (unit) out-of-sample performance. Table 2

Out-of-sample operational cost and its standard deviation with different choices of penalty coefficient γ and sparsity level d (bold-faced values represent the lowest cost under the same γ )

method Full model Exact method (d = 50) Exact method (d = 100) Exact method (d = 200) Exact method (d = 500) Exact method (d = 800)

γ = 0.5

γ=1

γ=5

γ = 10

mean

std

mean

std

mean

std

mean

std

4.5359 4.9280 4.7706 4.6522 4.5707 4.5374

0.1709 0.1842 0.1840 0.1748 0.1723 0.1696

4.0300 4.2086 4.0875 4.0460 4.0348 4.0281

0.1468 0.1672 0.1639 0.1507 0.1479 0.1475

3.7305 3.2440 3.3469 3.4718 3.6941 3.7296

0.1436 0.1154 0.1289 0.1363 0.1462 0.1435

3.7678 3.2073 3.3499 3.5306 3.7306 3.7670

0.1769 0.1246 0.1540 0.1586 0.1793 0.1757

We report the out-of-sample performance of the different methods in Table 2. Baseline approaches include the full model without variable selection, as well as our hard-cardinalityconstrained model (3) under various sparsity levels d. An ℓ2 regularization term with different coefficients γ is added into all models to mitigate overfitting. The numerical results

indicate that as the sparsity budget d increases, the out-of-sample performance of our variable selection framework approaches that of the full model. Moreover, with an appropriate choice of the ℓ2 regularization parameter and sparsity budget, our framework not only yields interpretable decisions but also attains the lowest out-of-sample operational costs among all methods considered. Table 3

Top 20 features selected using the exact algorithm with d = 100 on the real-world data, where the entries with

x represent the interaction between two features

Rank γ = 0.5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

γ=1

γ=5

γ = 10

n_days n_days n_days n_days n_items n_items n_items n_items transactions cluster_5_x_oil cluster_5_x_oil cluster_5 store_44 store_44_x_oil store_44_x_oil store_44 cluster_5 cluster_5 cluster_5 store_44_x_oil store_44_x_oil store_44 store_44 cluster_5_x_oil cluster_5_x_oil transactions transactions transactions cluster_8 cluster_8 month_9_x_store_9 city_Cuenca cluster_8_x_oil cluster_8_x_oil store_8_x_oil store_3 type_C_x_oil store_3 store_8 store_8 type_C store_8 store_3 store_9_x_oil store_3 store_8_x_oil city_Cuenca month_8_x_store_9 store_3_x_oil store_3_x_oil store_9 month_9_x_store_9 store_9 store_9 type_C month_2_x_store_44 store_8 type_C_x_oil store_9_x_oil type_C store_9_x_oil store_9_x_oil type_C_x_oil store_8_x_oil store_8_x_oil type_C store_3_x_oil city_Cuenca_x_oil city_Cuenca holiday_x_store_44 month_8_x_store_9 month_8_x_store_44 holiday_x_store_44 store_41_x_oil city_Cuenca_x_oil store_48_x_oil city_Quito store_41 cluster_10_x_oil store_9

We report the top 20 selected features under different configurations of (γ, d) in Tables 3 and 4. The features are ranked according to their averaged absolute weight values in 20 independent trials. From these selected features, we draw the following operational insights: (I) The number of days per week on which customers purchase hardware products (n_days) and the weekly number of active SKUs (n_items) emerge as the most persistent features. A possible explanation is that customer demand is largely driven by weekend shopping habits and the breadth of product assortment available in stores. (II) The historical weekly transaction volume (transactions) consistently appears as a strong predictor, as it is highly correlated with future customer demand. (III) The location and type of the store (such as those represented by the store_3, 8, 9, 44, C, city_Cuenca) play a critical role in predicting hardware sales. (IV) The interaction terms between store 9 and the months of August and September, as well

22 Top 20 features selected using the exact algorithm with d = 50 on the real-world data

Table 4

Rank γ = 0.5

γ=1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

n_days n_days n_days n_items n_items n_items cluster_5_x_oil cluster_5_x_oil store_44_x_oil store_44_x_oil store_44_x_oil cluster_5_x_oil store_44 cluster_5 cluster_5 transactions store_44 store_44 cluster_5 transactions transactions type_C_x_oil store_8 city_Cuenca type_C store_8_x_oil store_3 cluster_8 city_Cuenca store_8 cluster_8_x_oil store_3 store_8_x_oil store_9 store_9_x_oil month_9_x_store_9 store_9_x_oil city_Cuenca_x_oil month_8_x_store_9 store_8 store_9 store_9_x_oil store_8_x_oil type_C month_2_x_store_44 store_3 store_3_x_oil store_3_x_oil store_3_x_oil month_9_x_store_9 store_9 holiday_x_store_44 type_C_x_oil month_8_x_store_44 city_Cuenca cluster_8 type_C_x_oil city_Cuenca_x_oil store_48_x_oil store_48_x_oil

n_days n_items transactions store_44 cluster_5 store_44_x_oil cluster_5_x_oil type_C_x_oil type_C cluster_8 cluster_8_x_oil store_3 store_3_x_oil store_9 store_9_x_oil store_8 store_8_x_oil city_Cuenca city_Cuenca_x_oil city_Quito

γ=5

γ = 10

as between store 44 and February and August, suggest that the effects of store location and time on operational outcomes are heterogeneous, with certain stores exhibiting pronounced seasonal demand spikes.

8. Conclusion In this paper, we studied variable selection for the feature-based newsvendor problem under a hard cardinality constraint on the number of selected covariates. We establish both computational algorithms and statistical performance guarantees of this framework. The numerical results illustrate the trade-off among computational effort, support-recovery accuracy, and out-of-sample operating cost. Several directions remain open for future research. First, it is promising to extend this variable selection framework for nonlinear policy classes, such as the kernel-based decision rules [36]. Second, it is an interesting direction to design more efficient algorithms and hyperparameter selection procedures [12, 11]. Finally, extending the statistical analysis to misspecified, heavy-tailed, dependent, or distributionally robust settings would broaden the applicability of our framework.

References [1] Atamturk A, Gomez A (2025) Rank-one convexification for sparse regression. JMLR 1–50. [2] Bajwa WU, Mixon DG (2012) Group model selection using marginal correlations: The good, the bad and the ugly. 2012 50th Annual Allerton Conference on Communication, Control, and Computing, 494–501. [3] Ban GY, Rudin C (2019) The big data newsvendor: Practical insights from machine learning. Operations Research 67(1):90–108. [4] Belloni A, Chernozhukov V (2011) ℓ1 -penalized quantile regression in high-dimensional sparse models. The Annals of Statistics 39(1):82–130. [5] Bertsimas D, Kallus N (2020) From predictive to prescriptive analytics. Management Science 66(3):1025–1044. [6] Bertsimas D, Koduri N (2022) Data-driven optimization: A reproducing kernel hilbert space approach. Operations Research 70(1):454–471. [7] Bertsimas D, Pauphilet J, Van Parys B (2021) Sparse classification: a scalable discrete optimization perspective. Machine Learning 110(11-12):3177–3209. [8] Blanchet J, Kang Y, Murthy K (2019) Robust Wasserstein profile inference and applications to machine learning. Journal of Applied Probability 56(3):830–857. [9] Chang J, Yang L, Zhang Y, Zhou W (2025) Feature-rich, data-private: A sparse learning framework for the highdimensional newsvendor. SSRN 1–98. [10] Civril A, Magdon-Ismail M (2009) On selecting a maximum volume sub-matrix of a matrix and related problems. Theoretical Computer Science 410(47-49):4801–4811. [11] Cory-Wright R, Gómez A (2023) Optimal cross-validation for sparse linear regression. arXiv preprint arXiv:2306.14851 . [12] Cory-Wright R, Gómez A (2025) Stability regularized cross-validation. arXiv preprint arXiv:2505.06927 . [13] Dai S (2023) Variable selection in convex quantile regression: ℓ1 -norm or ℓ0 -norm regularization? European Journal of Operational Research 305(1):338–355. [14] Donti P, Amos B, Kolter JZ (2017) Task-based end-to-end model learning in stochastic optimization. Advances in neural information processing systems 30. [15] Dudley RM (1999) Uniform central limit theorems, volume 23 (Cambridge university press Cambridge). [16] Fan J, Li R (2001) Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American Statistical Association 96(456):1348–1360. [17] Favorita C (2017) Corporación favorita grocery sales forecasting. Kaggle, https://www.kaggle.com/competitions/favoritagrocery-sales-forecasting/data. [18] Günlük O, Linderoth J (2011) Perspective reformulation and applications. Mixed integer nonlinear programming, 61–89 (Springer). [19] Huang J, Horowitz JL, Ma S (2008) Asymptotic properties of bridge estimators in sparse high-dimensional regression models. The Annals of Statistics 36(2):587–613. [20] Huber J, Müller S, Fleischmann M, Stuckenschmidt H (2019) A data-driven newsvendor problem: From data to decision. European Journal of Operational Research 278(3):904–915.

24 [21] Kallus N, Mao X (2023) Stochastic optimization forests. Management Science 69(4):1975–1994. [22] Knight K (1998) Limiting distributions for L1 regression estimators under general conditions. The Annals of Statistics 26(2):755–770, URL http://dx.doi.org/10.1214/aos/1028144858. [23] Koenker R, Bassett G (1978) Regression quantiles. Econometrica 46(1):33–50. [24] Li Y, Dey SS, Xie W (2024) On sparse canonical correlation analysis. Advances in Neural Information Processing Systems 37:10707–10734. [25] Li Y, Xie W (2024) Beyond symmetry: best submatrix selection for the sparse truncated svd: Y. li, w. xie. Mathematical Programming 208(1):1–50. [26] Li Y, Xie W (2025) Exact and approximation algorithms for sparse principal component analysis. INFORMS Journal on Computing 37(3):582–602. [27] Liyanage LH, Shanthikumar JG (2005) A practical inventory control policy using operational statistics. Operations research letters 33(4):341–348. [28] Notz PM, Pibernik R (2022) Prescriptive analytics for flexible capacity management. Management Science 68(3):1756– 1775. [29] Oroojlooyjadid A, Snyder LV, Takáč M (2020) Applying deep learning to the newsvendor problem. IISE Transactions 52(4):444–463. [30] Qi M, Grigas P, Shen ZJ (2025) Integrated conditional estimation-optimization. Operations Research . [31] Sadana U, Chenreddy A, Delage E, Forel A, Frejinger E, Vidal T (2025) A survey of contextual optimization methods for decision-making under uncertainty. European Journal of Operational Research 320(2):271–289. [32] Serrano B, Minner S, Schiffer M, Vidal T (2024) Bilevel optimization for feature selection in the data-driven newsvendor problem. European Journal of Operational Research 315(2):703–714. [33] Tibshirani R (1996) Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology 58(1):267–288. [34] Tropp JA (2010) User-friendly tail bounds for sums of random matrices. Foundations of Computational Mathematics 12:389–434. [35] Wainwright MJ (2019) High-dimensional statistics: A non-asymptotic viewpoint, volume 48 (Cambridge university press). [36] Wang J, Dey SS, Xie Y (2023) Variable selection for kernel two-sample tests. arXiv preprint arXiv:2302.07415 . [37] Xie H, Huang J (2009) SCAD-penalized regression in high-dimensional partially linear models. The Annals of Statistics 37(2):673–696. [38] Xie W, Deng X (2020) Scalable algorithms for the sparse ridge regression. SIAM Journal on Optimization 30(4):3359– 3386. [39] Zhang CH (2010) Nearly unbiased variable selection under minimax concave penalty. The Annals of Statistics 38(2):894–942. [40] Zhang Y, Gao J (2017) Assessing the performance of deep learning algorithms for newsvendor problem. International conference on neural information processing, 912–921.

ec1

Appendix for “Variable Selection for Feature-Based Newsvendor” EC.1. Proof of Proposition 1 We prove the NP-hardness result by making the reduction from Exact Cover by

3-Sets (X3C) to Problem (3). The formal description of (X3C) is as follows. Let Q be a ground set consisting of m := 3d elements, and C = {C1 , . . . , CD } be the collection of subsets of Q such that |Ci | = 3 for i ∈ [D]. Problem (X3C) aims to determine whether there exists a subset of C , denoted as C ′ such that every element in Q occurs exactly once in C ′ . Denote by A ∈ Rm×D the incidence matrix of Problem (X3C) such that A[i, j] = 1 if the i-th element of Q belongs to the subset Cj and otherwise A[i, j] = 0, ∀i, j. We define a scalar M0 m . We construct a special instance of Problem (3) by taking n = 2m and such that M0 > 3γ(b∧h)

the i-th observation as

(xi , yi ) =

 

(M0 A[i, :]⊤ , M0 ) if 1 ≤ i ≤ m,

(−M A[i − m, :]⊤ , −M ) if m + 1 ≤ i ≤ n. 0 0

For weight vector ω ∈ RD , define the error vector r = 1 − Aω ∈ Rm . For any scalar x, define

x+ = max{x, 0} and x− = max{−x, 0}. Define b ∧ h = min(b, h). Thus, the fitting error part can be bounded as n i 1 Xh b · max{0, yi − ω ⊤ xi − c} + h · max{0, ω ⊤ xi + c − yi } n i=1 m i 1 Xh b · (M0 r[i] − c)+ + h · (M0 r[i] − c)− + b · (−M0 r[i] − c)+ + h · (−M0 r[i] − c)− = 2m i=1 m i b∧h Xh ≥ |M0 r[i] − c| + |M0 r[i] + c| 2m i=1

≥

(b ∧ h)M0 ∥r∥1 , m

where the last inequality is based on the relation |u + v| + |u − v| ≥ 2|u|. Let (ω, c) be a feasible solution to Problem (3) with ∥ω∥0 ≤ d. We first show that F (ω, c) ≥

d/(2γ). • Assume that ∥ω∥22 ≤ d, then ⊤

|1 ω| =

X i: ω[i]̸=0

ω[i] ≤

s X i: ω[i]̸=0

s X √ 1 ω 2 [i] ≤ d∥ω∥2 , i: ω[i]̸=0

(EC.1)

ec2

where the first inequality uses Cauchy-Schwarz inequality. It follows that

|1⊤ ω|2 = (d − 1⊤ ω) · 1 + 1⊤ ω/d d ≤ 2(d − 1⊤ ω) 2 2 = 1⊤ r ≤ ∥r∥1 , 3 3

d − ∥ω∥22 ≤ d −

where the first inequality is by (EC.1), the second inequality is because |1⊤ ω| ≤ √ d∥ω∥2 ≤ d, and the last equality is because each column of A has exactly three ones and 1⊤ r = m − 1⊤ Aω = m − (A⊤ 1)⊤ ω = 3d − 3(1⊤ ω).

(EC.2)

2 d d − ∥r∥1 ≥ , 3 2γ

(EC.3)

then, the objective

1 (b ∧ h)M0 ∥r∥1 + F (ω, c) ≥ m 2γ where the last inequality is by the choice of M0 . • Otherwise, for ∥ω∥22 > d, the objective value

F (ω, c) ≥

1 d ∥ω∥22 > . 2γ 2γ

We now show that Problem (X3C) can be reduced to Problem (3). Assume we solve Probd and denote its optimal solution as (ω, c). lem (3) with optimal value upper bounded by 2γ

This happens only if ∥ω∥22 ≤ d, and the inequalities in (EC.3) all hold with equalities. It implies that r = 0 = 1 − Aω . By (EC.2), it holds that 1⊤ ω = d. By (EC.1), it holds that ∥ω∥22 = d. Then, all inequalities in (EC.1) always hold with equalities, and therefore ω has exactly d nonzero entries. The exact cover of (X3C) is given by {Cj : ω[j] ̸= 0}. Conversely, suppose Problem (X3C) has an exact cover. Let ω ∈ RD be the associated weight vector such that ω[i] = 1 if the i-th element of C is selected to form the cover and otherwise

ω[i] = 0. Then the incidence matrix A satisfies that r := 1 − Aω = 0 and ω satisfies ∥ω∥0 = d and ∥ω∥22 = d. The objective value of the feasible solution (ω, 0) to Problem (3) equals m i d 1 Xh F (ω, 0) = b · (M0 r[i])+ + h · (M0 r[i])− + b · (−M0 r[i])+ + h · (−M0 r[i])− + 2m i=1 2γ

d . 2γ

It justifies that the optimal value of Problem (3) is bounded by d/(2γ). The proof is completed.

ec3

EC.2. Proof of Theorem 1 Both Big-M and perspective reformulation can be written as the following form of optimization:

min

ω,c,s,µ

s.t.

n i 1 1 Xh b max{0, yi − ω ⊤ xi − c} + h max{0, ω ⊤ xi + c − yi } + ∥µ∥1 n i=1 2γ

(EC.4)

(ω, c, s, µ) ∈ V Big-M or V Pers , X s[j] ≤ d, j∈[D]

where

    2 ω[i] ≤ µ[i], i ∈ [D] , × R : V Big-M = (ω, s, µ, c) ∈ RD × {0, 1}D × RD +  |ω[i]| ≤ M [i] s[i], i ∈ [D]      2 ω[i] ≤ µ[i]s[i], i ∈ [D] Pers D D D V = (ω, s, µ, c) ∈ R × {0, 1} × R+ × R : .  |ω[i]| ≤ M [i] s[i], i ∈ [D]  These two problems are equivalent because V Big-M = V Pers . Indeed, when the binary variable

s[i] = 0, it holds that ω[i] = 0 and the first constraints appeared in V Big-M and V Pers are redundant. Otherwise, the first constraints in V Big-M and V Pers are equivalent. For Part (II), if considering the continuous relaxation of these two problems, the binary constraints s[i] ∈ {0, 1} in V Big-M and V Pers are replaced by s[i] ∈ [0, 1] for i ∈ [D]. Denoted the continuous relaxation of these two sets as V

Big-M

and V

Pers

. It holds that V

Pers

⊆V

Big-M

ω[i]2 ≤ µ[i]s[i] implies ω[i]2 ≤ µ[i]. Therefore, the optimal value of continuous relaxation of perspective reformulation is larger than that of Big-M reformulation.

EC.3. Proof of Theorem 2 For Part (I), we bound the support size of solution in Algorithm 1. Let (b ω, b c, sb, µ b) be the optimal solution to the continuous relaxation of (Perspective). Let (ω̃, c̃, s̃) be the output of Algorithm 1. The support size of ω̃ is bounded by the sum of random variables N := P s[j]). Its expected value j∈[D] s̃[j], where s̃[j] ∼ Bernoulli(b

E [N ] =

X j∈[D]

sb[j] ≤ d.

ec4

For t > 0, we compute the moment-generating function of N :

MN (t) = E[etN ] =

E[ets̃[j] ]

j∈[D]

Y 1 + (b s[j]et − sb[j]) ≤ exp sb[j]et − sb[j]

j∈[D]

 = exp (et − 1)

 X

sb[j] ≤ exp(d(et − 1)),

j∈[D]

where the first inequality is based on the relation that 1 + x ≤ ex , ∀x, and the last inequality P is because j∈[D] sb[j] ≤ d. For ∆ > 0, applying the Markov inequality gives

MN (t) P (N ≥ (1 + ∆)d) = P etN ≥ et(1+∆)d ≤ t(1+∆)d e ≤ exp(d(et − 1) − t(1 + ∆)d). We choose t such that the upper bound above is as sharp as possible, yielding the optimal choice t = log(1 + ∆). Then, it holds that

P (N ≥ (1 + ∆)d) ≤ exp(−d((1 + ∆) log(1 + ∆) − ∆)). 2∆ For ∆ ≥ 0, it is easy to verify that log(1 + ∆) ≥ 1+∆ , and then d∆2 P (N ≥ (1 + ∆)d) ≤ exp − . 2+∆ q 2 Taking ∆ = 2ξ + 2ξ + ξ4 with ξ := log(2/δ)/d yields

δ (EC.5) P (∥ω̃∥0 ≥ (1 + ∆)d) ≤ P (N ≥ (1 + ∆)d) ≤ . 2 p We replace this bound with a conservative one (∆ = log(2/δ)/d + 2 log(2/δ)/d) to simplify the notation in the theorem statement. For Part (II), we first reformulate the continuous relaxation of Problem (Perspective):

min

ω,c,s,µ

s.t.

n i 1 1 Xh ⊤ ⊤ b max{0, yi − ω xi − c} + h max{0, ω xi + c − yi } + ∥µ∥1 n i=1 2γ

(ω, c, s, µ) ∈ RD × R × [0, 1]D × RD +, X s[j] ≤ d, j∈[D]

ω[i]2 ≤ µ[i]s[i], i ∈ [D].

(EC.6)

ec5 2

The variable µ above can be eliminated by taking µ[i] = ω[i] (by default, we take µ[i] = 0 if s[i]

ω[i] = s[i] = 0), yielding the reformulation min ω,c,s

s.t.

n D i X 1 Xh ω[i]2 ⊤ ⊤ b max{0, yi − ω xi − c} + h max{0, ω xi + c − yi } + n i=1 2γs[i] i=1

(EC.7)

(ω, c, s) ∈ RD × R × [0, 1]D , X s[j] ≤ d. j∈[D]

By fixing s and taking the duality of the remaining optimization problem, we obtain the min-max reformulation of the continuous relaxation: γ 1 ⊤ ⊤ ⊤ XDiag(s)X z − z y . min g(s) = maxP −z 2 D, z[i]∈[−b,h],i∈[n], n z[i]=0 2n n s∈[0,1] i =1 P

(EC.8)

j∈[D ] s[j]≤d

It is also clear that optval(3) = min {g(s) : s ∈ {0, 1}D ,

j∈[D] s[j] ≤ d}. By definition,

F (ω̃, c̃) − optval(3) ≤ F (ω̃, c̃) − g(b s) = g(e s) − g(b s),

(EC.9)

where the first inequality is because sb is the optimal solution obtained from continuous relaxation, and the second inequality is by strong duality. Next, we provide the upper bound on g(s̃) − g(b s). Define the inner objective of g(s) as

g(s, z) := −z ⊤

1 ⊤ ⊤ XDiag(s)X z − z y. 2n2 n

Let z̃ be the maximizer of the inner problem of g(s̃), then

g(s̃) − g(b s) ≤ g(s̃, z̃) − g(b s, z̃) γ = − 2 z̃ ⊤ XDiag(s̃ − sb)X ⊤ z̃ 2n γ 2 ⊤ ≤ 2 ∥z̃∥2 λmax XDiag(b s − s̃)X 2n γ(b ∨ h)2 ≤ λmax XDiag(b s − s̃)X ⊤ , 2n

(EC.10)

where the first inequality is because z̃ is a feasible solution (but perhaps suboptimal) to the inner problem of g(b s), the second inequality is based on the Cauchy-Schwarz inequality, and the last one is because ∥z̃∥22 ≤ n(b ∨ h)2 . It suffices to build the concentration bound P ⊤ on λmax XDiag(b s − s̃)X ⊤ = λmax R s[i] − s̃[i])X:,i X:,i . Following the i∈[D] i with Ri = (b similar argument as in [38, Proof of Lemma 3], it can be shown that E[Ri ] = 0, λmax (Ri ) ≤ θ1

ec6

and ∥

2 i∈[D] E[Ri ]∥2 ≤ θ1 θd . Applying the matrix concentration bound in [34, Theorem 1.4]

yields

n o ⊤ P λmax XDiag(b s − s̃)X ≥ t ≤ n exp − Taking t = 32 θ1 log(2n/δ) +

t2 . 2θ1 θd + 2/3θ1 t

(EC.11)

2θ1 θd log(2n/δ) ensures the right-hand-side of (EC.11) being

bounded by δ/2. Combining the inequalities in (EC.9), (EC.10), and (EC.11) yields p γ(b ∨ h)2 2 P F (ω̃, c̃) − optval(3) ≤ · θ1 log(2n/δ) + 2θ1 θd log(2n/δ) 2n 3 2 p ⊤ ≥P λmax XDiag(b s − s̃)X ≤ θ1 log(2n/δ) + 2θ1 θd log(2n/δ) 3 δ ≥1 − . 2 The choice of γ in the main statement of Theorem 2 further ensures that

P {F (ω̃, c̃) − optval(3) ≤ ϵ} ≥ 1 − δ/2.

(EC.12)

Applying the union bound of the concentration results in (EC.5) and (EC.12) gives the desired result.

EC.4. Proof of Theorem 4 Let us define the empirical cost

Rn (ω, c) =

n i 1 Xh b · max{0, yi − ω ⊤ xi − c} + h · max{0, ω ⊤ xi + c − yi } . n i=1

and the population cost

R(ω, c) = E(X,Y ) b · max{0, Y − ω ⊤ X − c} + h · max{0, ω ⊤ X + c − Y } . We rely on the following two technical lemmas that are helpful to show the proof of Theorem 4. Lemma EC.1 (Global Risk Lower Bound). Under Assumption 1, for any ∆ := (u⊤ , t)⊤ ∈ RD+1 with ∥u∥0 ≤ 2d, it holds that

R(ω ∗ + u, c∗ + t) − R(ω ∗ , c∗ ) ≥

(b + h)f κr0 · ∥∆∥22 . 2(r0 + K∥∆∥2 )

ec7

Proof. By definition,

R(ω ∗ + u, c∗ + t) − R(ω ∗ , c∗ ) =(b + h)E(X,Y ) ρτ (Y − (ω ∗ + u)⊤ X − (c∗ + t)) − ρτ (Y − (ω ∗ )⊤ X − c∗ ) h i e − ρτ (ε) , =(b + h)E(ε,X) ρτ (ε − ∆⊤ X) where the first equality is by the definition of quantile function ρτ (u) = u(τ − 1{u<0} ), and the e = (X ⊤ , 1)⊤ , ∆ = second one is based on the model assumption (7) and the expressions that X

(u⊤ , t)⊤ . Based on the Knight’s identity [22], it holds that h i e − ρτ (ε) E(ε,X) ρτ (ε − ∆⊤ X) " # Z ∆⊤ Xe e − 1{ε<0} ) + 1{ε≤s} − 1{ε<0} ds =E(ε,X) −∆⊤ X(τ

(EC.13)

"Z =E(ε,X)

e ∆⊤ X

1{ε≤s} − 1{ε<0} ds = EX

e ∆⊤ X

Fε (s) − Fε (0) ds ,

where the second and last equalities are due to the assumption that ε and X are independent together with the definition of the cumulative distribution function (cdf) Fε (s) := P(ε ≤ s). R ∆⊤ Xe e ≤ r0 , by Fε (s) − Fε (0) ds as follows. When |∆⊤ X| We provide lower bound on 0

Assumption 1(III), Z ∆⊤ Xe

Z |∆⊤ X| e e 2 f (∆⊤ X) Fε (s) − Fε (0) ds ≥ f s ds = , 2 0 0 Rs Rs where the inequality is because Fε (s) − Fε (0) = 0 fε (t) dt ≥ 0 f dt = f s for s ≥ 0. When e > r0 > 0, it holds that ∆⊤ X Z ∆⊤ Xe

Z r0 Z ∆⊤ Xe Fε (s) − Fε (0) ds ≥ Fε (s) − Fε (0) ds + Fε (r0 ) − Fε (0) ds

≥

f r02 2

e − r0 )f r0 ≥ + (∆⊤ X

0 ⊤ e

f r0 ∆ X , 2

e into two parts and by the where the first inequality is by separating the integral on [0, ∆⊤ X] monotonicity of the cdf, the second inequality is by Assumption 1(III), and the last one is R ⊤ e > r0 . It can be shown using the similar argument that ∆ Xe Fε (s) − by the relation ∆⊤ X 0 ⊤ e ⊤ e Fε (0) ds ≥ −f r0 ∆ X/2 for ∆ X < −r0 < 0. Combining these cases yields the global lower bound

Z ∆⊤ Xe 0

e 2 f r0 (∆⊤ X) e e ≥ Fε (s) − Fε (0) ds ≥ f |∆⊤ X|/2 · min{r0 , |∆⊤ X|} . e 2(r0 + |∆⊤ X|)

ec8

Substituting this lower bound to (EC.13) yields

R(ω ∗ + u, c∗ + t) − R(ω ∗ , c∗ ) " # " # e 2 e 2 (b + h)f r0 (∆⊤ X) (b + h)f r0 κ∥∆∥22 (b + h)f r0 (∆⊤ X) ≥EX ≥ EX ≥ , e 2(r0 + ∥∆∥2 K) 2(r0 + ∥∆∥2 K) 2(r0 + |∆⊤ X|) where the second and last inequalities leverage Assumptions 1(IV) and 1(V), respectively.

□

Lemma EC.2 (Uniform Concentration Bound on Risk). Under Assumption 1, there exists a constant C0 depending on b, h, K such that, with probability at least 1 − δ , simultaneously for every ∆ := (u⊤ , t)⊤ ∈ RD+1 with ∥u∥0 ≤ 2d, it holds that Rn (ω ∗ + u, c∗ + t) − Rn (ω ∗ , c∗ ) − R(ω ∗ + u, c∗ + t) − R(ω ∗ , c∗ ) r d log(eD/d) + 1 + log(2/δ) . ≤ C0 ∥∆∥2 · n Proof. The inequality for ∆ = 0 holds trivially. We focus on the case where ∆ ̸= 0 instead. Define the function ℓ(y, z) = b max{0, y − z} + h max{0, z − y}, and the function class ℓ(y, x⊤ (ω ∗ + u) + c∗ + t) − ℓ(y, x⊤ ω ∗ + c∗ ) ⊤ ⊤ H = (x, y) 7→ : ∆ ̸= 0, ∆ := (u , t) , ∥u∥0 ≤ 2d . ∥∆∥2 It suffices to provide concentration bound on n

1X sup EX,Y [h(X, Y )] − h(xi , yi ) , n i=1 h∈H which can be shown based on the covering number argument. Each element in the function class H can be indexed with two terms: r := ∥∆∥2 and v := ∆ . Hence, we reformulate ∥∆∥2

  ⊤ ∗ ∗ ⊤ ∗ ∗ ℓ(y, x (ω + rv ) + c + rv ) − ℓ(y, x ω + c )   u t (x, y) 7→ hr,v (x, y) := , r H= .   r > 0, v = (v ⊤ , v )⊤ ∈ RD+1 , ∥v∥ = 1, ∥v ∥ ≤ 2d  t 2 u 0 u Define the empirical Rademacher complexity " # n X 1 b n (H) := Eσ ,i∈[n] sup R σi h(xi , yi ) (xi , yi ), i ∈ [n] , i h∈H n i=1

b n (H) where σi , i ∈ [n] are i.i.d. Rademacher random variables. We provide upper bound on R using the following procedure.

ec9

• For fixed v ∈ RD+1 , denote the subset of H as n o Hv = hr,v : r > 0 . Define the empirical L2 pseudo-metric as ( n ) 2 1/2 1 X dn (f, g) = f (xi , yi ) − g(xi , yi ) . n i=1 Let r1 , . . . , rM with 0 < r1 < · · · < rM index an δ -packing of Hv under the metric dn (·). Therefore, n M −1

M −1 X

2 1XX (M − 1)δ < hrj ,v (xi , yi ) − hrj +1 ,v (xi , yi ) . dn (hrj ,v , hrj +1 ,v ) = n i=1 j=1 j=1 2

According to the formulation of hr,v and by the convexity of ℓ(y, ·), it is easy to show that

hr,v (x, y) is non-decreasing in r, and therefore n M −1 2 1XX hrj ,v (xi , yi ) − hrj +1 ,v (xi , yi ) (M − 1)δ < n i=1 j=1 "M −1 # n 2 1X X hrj +1 ,v (xi , yi ) − hrj ,v (xi , yi ) ≤ n i=1 j=1 2

2 1 X = hrM ,v (xi , yi ) − hr1 ,v (xi , yi ) ≤ 4(b ∨ h)2 K 2 , n i=1 where the last inequality uses the fact that for any (x, y) and r > 0,

|hr,v (x, y)| ≤ (b ∨ h) ·

|x⊤ (rvu ) + rvt | ≤ (b ∨ h)∥(x⊤ , 1)⊤ ∥2 ∥v∥2 ≤ (b ∨ h)K. r

The packing number M ≤ 1 + 4(b ∨ h)2 K 2 /δ 2 . Therefore, the covering number of Hv under dn (·) is bounded:

N δ, Hv , dn ≤ 1 + 4(b ∨ h)2 K 2 /δ 2 . • Next, we study the covering number of the set n o V = v = (vu⊤ , vt )⊤ ∈ RD+1 : ∥v∥2 = 1, ∥vu ∥0 ≤ 2d .

n For fixed support S ⊆ [D] with |S| ≤ 2d, the δ -covering number of VJ = v = (vu⊤ , vt )⊤ ∈ o RD+1 : ∥v∥2 = 1, supp(vu ) = J under ∥ · ∥2 is at most (1 + 2/δ)2d+1 . By taking the union of all sparse supports, the covering number   min{2d,D} X D  N (δ, V, ∥ · ∥2 ) ≤  (1 + 2/δ)2d+1 . k k=0

ec10

• For fixed r > 0 and any two unit vectors v = (vu , vt ), v ′ = (vu′ , vt′ ), it holds that

hr,v (x, y) − hr,v′ (x, y) ≤ (b ∨ h) ·

|x⊤ (rvu − rvu′ ) + rvt − rvt′ | ≤ (b ∨ h)K∥v − v ′ ∥2 . r

For any hr,v ∈ H, we select v ′ in the δ -cover of V and hr′ ,v′ in the δ -cover of Hv′ . It follows that

# n 2 1/2 1 X dn (hr,v , hr′ ,v′ ) = hr,v (xi , yi ) − hr′ ,v′ (xi , yi ) n i=1 # " n n 2 2 X 2 1/2 2 X hr,v (xi , yi ) − hr,v′ (xi , yi ) + hr′ ,v′ (xi , yi ) − hr,v′ (xi , yi ) ≤ n i=1 n i=1 1/2 = (2 + 2(b ∨ h)2 K 2 )1/2 δ. ≤ 2(b ∨ h)2 K 2 δ 2 + 2δ 2 "

Therefore,

N (2 + 2(b ∨ h)2 K 2 )1/2 δ, H, dn ≤ N δ, Hv , dn · N (δ, V, ∥ · ∥2 ). It further implies that

log N δ, H, dn

  min{2d,D} 2 2 2 2 X 8(b ∨ h) K (1 + (b ∨ h) K ) D  ≤ log 1 + + log  2 δ k k=0 " # √ 2 2(1 + (b ∨ h)2 K 2 )1/2 + (2d + 1) log 1 + δ eD 1 ≤C1 · d log + 1 + d log , d δ

where C1 is a constant depending on b, h, K . As each element in H is uniformly bounded by (b ∨ h)K , by Dudley’s entropy integral bound [15], the Rademacher complexity r Z C2 (b∨h)K Rn (H) ≤ √ log N δ, H, dn dδ n 0 s d log eD +1 d ≤ C3 , n where C2 , C3 are some constants depending on b, h, K . Based on the symmetrization and McDiarmid’s inequalities, it holds with probability at least

1 − δ that s

1X h(xi , yi ) − E[h(X, Y )] ≤ C4 sup h∈H n i=1

eD d

n r

≤ C0

d log

r + (b ∨ h)K

d log(eD/d) + 1 + log(2/δ) , n

2 log(2/δ) n

ec11

where C0 , C4 are some constants depending on b, h, K . Re-arranging the inequality above yields the desired result.

□

Finally, we provide the bound on the residual error   ∗ ω bn − ω  an :=  b cn − c∗

. 2

Proof of Theorem 4. As ∥b ωn − ω ∗ ∥0 ≤ 2d, by Lemma EC.1, we obtain the lower bound

R(b ωn , b cn ) − R(ω ∗ , c∗ ) ≥

(b + h)f κr0 2 ·a . 2(r0 + Kan ) n

(EC.14)

Next, we derive its upper bound. By the definition of (b ωn , b cn ), it holds that

Rn (b ωn , b cn ) +

1 1 ∥b ωn ∥22 ≤ Rn (ω ∗ , c∗ ) + ∥ω ∗ ∥22 . 2γn 2γn

Leveraging this relation and re-arranging the terms yields

R(b ωn , b cn ) − R(ω ∗ , c∗ ) h i h i ∗ ∗ ∗ ∗ ∗ ∗ = Rn (b ωn , b cn ) − Rn (ω , c ) − Rn (b ωn , b cn ) − Rn (ω , c ) − R(b ωn , b cn ) − R(ω , c ) h i h i 1 ∗ 2 2 ∗ ∗ ∗ ∗ ∥ω ∥2 − ∥b ωn ∥2 − Rn (b ωn , b cn ) − Rn (ω , c ) − R(b ωn , b cn ) − R(ω , c ) . ≤ 2γn The first term above can be further bounded as follows:

i ⟨ω ∗ , ω ∗ − ω bn ⟩ ∥ω ∗ ∥2 an 1 h ∗ 2 ∥ω ∥2 − ∥b ωn ∥22 ≤ ≤ . 2γn γn γn By Lemma EC.2, with probability at least 1 − δ , the second term above can be bounded as r d log(eD/d) + 1 + log(2/δ) C0 an . n Combining these two relations implies the following relation holds with probability at least

1 − δ: R(b ωn , b cn ) − R(ω ∗ , c∗ ) ≤ an ·

∥ω ∗ ∥2 + C0 γn

d log(eD/d) + 1 + log(2/δ) n

(EC.15)

Comparing this with (EC.14) yields

(b + h)f κr0 2 · a ≤ an · 2(r0 + Kan ) n

∥ω ∗ ∥2 + C0 γn

d log(eD/d) + 1 + log(2/δ) n

ec12

Specify 2γ1n ≲

d log(eD/d) , then there exists constant C1 depending on C0 , b, h, K, ∥ω ∗ ∥2 such n

that

(b + h)f κr0 · an ≤ C1 2(r0 + Kan )

d log(eD/d) + 1 + log(2/δ) . n

(EC.16)

Assume n is sufficiently large such that r d log(eD/d) + 1 + log(2/δ) (b + h)f κr0 ≤ , C1 n 4K then re-arranging (EC.16) yields q d log(eD/d)+1+log(2/δ) C1 r0 n q an ≤ (b+h)f κr0 d log(eD/d)+1+log(2/δ) − C1 K 2 n q d log(eD/d)+1+log(2/δ) r C1 r0 n 4C1 d log(eD/d) + 1 + log(2/δ) · ≤ = (b+h)f κr0 (b + h)f κ n 4

holds with probability at least 1 − δ .

□

EC.5. Proof of Corollary 1 By definition,

min

∥ω∥0 ≤d,c

R(ω, c) = min

∥ω∥0 ≤d,c

(b + h)E(X,Y ) ρτ (Y − (ω ⊤ X + c)) .

By the likelihood model, for each x, (ω ∗ )⊤ x + c∗ is the τ -quantile of Y | X = x, and therefore the optimization problem above takes the optimal solution (ω ∗ , c∗ ). Recall from the proof of Theorem 4 that for sufficiently large n, with probability at least 1 − δ , the following two relations hold:

∗

d log(eD/d) + 1 + log(2/δ) , n r 4C1 d log(eD/d) + 1 + log(2/δ) an ≤ · . (b + h)f κ n ∗

R(b ωn , b cn ) − R(ω , c ) ≤ an C1 ·

Combining these two relations yields

R(b ωn , b cn ) − R(ω ∗ , c∗ ) ≤

4C12 d log(eD/d) + 1 + log(2/δ) · . (b + h)f κ n

ec13

EC.6. Proof of Corollary 2 Under the event described in Corollary 2, for each j ∈ S ∗ , it holds that   ∗ ω bn − ω  > 0. |b ωn [j]| ≥ |ω ∗ [j]| − |b ωn [j] − ω ∗ [j]| ≥ |ωj∗ | −  b cn − c∗ 2

Therefore, S ∗ ⊆ supp(b ωn ). On the other hand, by the feasibility of ω bn , it holds that

|supp(b ωn )| ≤ d = |S ∗ |. This fact, together with S ∗ ⊆ supp(b ωn ) implies that S ∗ = supp(b ωn ).

Variable Selection for Feature-Based Newsvendor

Analyze & explore related research in Conceptio

Related documents