Sparsity Regularized and Robust Mean Variance Portfolio Selection Under Ellipsoidal Uncertainty
arXiv:2609.11749v1 [math.OC] 10 Sep 2026
Deniz Akkaya∗
Emre Can Yayla∗
Buse Şen†
Mustafa Ç. Pınar∗
Abstract We investigate mean–variance portfolio selection with an ℓ0 -penalty to promote sparsity in asset allocations. Uncertainty in the mean return vector is incorporated through an ellipsoidal uncertainty set, yielding a robust sparse optimization framework. We characterize the structure of both local and global minimizers and exploit these properties in the risk minimization and return maximization formulations. Building on this structural insight, we develop a branchand-bound algorithm tailored to the resulting robust sparse portfolio problems, together with a new pruning rule that can discard exponentially many candidate portfolios in a single step. Extensive computational experiments on real market data, together with comparisons against a mixed-integer second-order cone programming solver, demonstrate the effectiveness and competitiveness of the proposed approach. Keywords: Mean-variance portfolio, robust optimization, regularization, sparsity, ℓ0 -norm, Branch-and-Bound Mathematics Subject Classification (MSC 2020): 90C11, 91G10, 90C17, 90C26, 90C57
1
Introduction
Mean-variance portfolio selection, originally introduced by [19], forms the cornerstone of modern portfolio theory. In this framework, an investor determines portfolio weights by balancing expected return against risk measured through variance, thereby tracing the efficient frontier. Despite its conceptual elegance and widespread adoption, the classical mean-variance model is known to be highly sensitive to estimation errors in the input parameters, particularly in the expected return vector. As emphasized in the literature, small perturbations in estimated returns may lead to large and economically unintuitive changes in optimal portfolios, a phenomenon often attributed to the amplification of estimation noise [20]. To mitigate this instability, robust optimization has emerged as a systematic approach for incorporating parameter uncertainty directly into optimization models [5, 6, 14]. In robust portfolio optimization, uncertain parameters such as expected returns are assumed to lie within prescribed uncertainty sets, and portfolio decisions are made to perform well under the worst-case realization within these sets. A prominent modeling choice is the use of ellipsoidal uncertainty sets for the mean return vector, which lead to tractable reformulations and admit appealing statistical interpretations. In particular, [16] demonstrated that robust mean-variance portfolio problems with ellipsoidal uncertainty can be reformulated as convex optimization problems and that the resulting portfolios exhibit improved stability and out-of-sample performance. Subsequent contributions, including [8], explored alternative uncertainty structures and their computational implications. In ∗ †
Department of Industrial Engineering, Bilkent University, Ankara, Türkiye Risk Analytics and Optimization Chair, EPFL, Lausanne, Switzerland
1
addition, distributionally robust optimization has provided a broader framework for modeling uncertainty in financial decision making, where ambiguity in probability distributions is explicitly incorporated [13, 15, 25]. Among various robust formulations, we follow the analysis of robust portfolio models developed in [23] as a foundation for our study. While robustness addresses estimation risk, practical portfolio construction typically involves additional structural considerations. In many real-world applications, investors restrict the number of assets held in a portfolio due to transaction costs, liquidity constraints, monitoring costs, or regulatory requirements. Such considerations naturally lead to sparse portfolio models, where the number of nonzero positions is explicitly controlled. Early work on cardinality-constrained portfolio optimization highlighted the computational challenges arising from these restrictions [9, 12, 18]. More recently, sparsity-inducing formulations and scalable algorithms have attracted considerable attention in the literature. Contributions such as [7, 10, 17, 26] demonstrate the practical and computational advantages of sparse portfolios. A direct way to enforce sparsity is through cardinality constraints or ℓ0 -regularization, where the ℓ0 -term counts the number of selected assets. However, the inclusion of such a term leads to nonconvex and combinatorial optimization problems that are challenging to solve. In [22], a comprehensive analysis of sparse portfolio optimization is presented, including structural properties of optimal portfolios and an efficient enumeration-based branching algorithm. In addition, recent studies continue to explore models that combine sparsity and robustness in portfolio construction [1, 27]. Despite the extensive literature on robust portfolio optimization and the growing body of work on sparse portfolio selection, the integration of ellipsoidal mean uncertainty with exact ℓ0 regularization remains relatively unexplored. Most robust formulations focus primarily on mitigating estimation risk but do not explicitly control portfolio cardinality. Conversely, sparse portfolio models often rely on deterministic parameter estimates and do not explicitly account for uncertainty in expected returns. A unified treatment that simultaneously addresses estimation risk and enforces exact sparsity leads to a challenging class of robust mixed-integer quadratic optimization problems. Motivating this integration, preliminary out-of-sample tests on real equity market data show that adding ellipsoidal robustness to a sparse mean-variance model can improve risk-adjusted performance over the purely sparse benchmark. In this paper, we study sparse mean-variance portfolio selection under ellipsoidal uncertainty in the mean return vector. We consider both risk minimization and return maximization variants and incorporate an ℓ0 -penalty to promote portfolios with a prescribed level of sparsity. The ellipsoidal uncertainty set is consistent with the robust optimization framework developed in [6, 16, 23], while the ℓ0 -regularization follows the exact sparsity-inducing perspective studied in [2, 3, 21] and recent sparse portfolio optimization research presented in [22]. Our contributions are both theoretical and computational. On the theoretical side, we characterize the structure of local and global minimizers of the robust sparsity-penalized problem. We analyze how the interaction between the quadratic risk term, the worst-case mean adjustment induced by the ellipsoidal uncertainty set, and the discontinuous ℓ0 -penalty determines the structure of optimal portfolios. These results extend structural analyses known for deterministic sparse mean-variance models to a robust setting and clarify the role of the robustness parameter in shaping portfolio composition. On the computational side, we develop a tailored branch-and-bound framework that exploits problem-specific lower and upper bounds derived from the analytical structure of the model. These bounds enable effective pruning of candidate supports and provide a solid basis for warm-start heuristics. We conduct extensive computational experiments on real financial data and benchmark our approach against mixed-integer second-order conic formulations. The numerical results show that the proposed method is computationally effective, especially on larger instances, and often 2
produces high-quality sparse robust portfolios with substantially reduced running times. Overall, the paper contributes both structural analysis and algorithmic developments for robust sparse mean-variance portfolio optimization. The main contributions of this paper are as follows. • We propose a robust sparse mean-variance portfolio model that integrates ellipsoidal uncertainty in expected returns with exact ℓ0 -regularization, providing a unified treatment of robustness and sparsity. • We develop a structural analysis of the resulting nonconvex problem, including support-wise decomposition, identification of local minimizers and results on existence and properties of global minimizers. • We establish explicit bounds on optimal solutions and show that sparsity levels can be controlled via the regularization parameter through thresholding results. • We design a tailored branch-and-bound algorithm that turns these structural insights into bounding and warm-start strategies, and show on real data that it is computationally effective. A key ingredient is an additional pruning step that refines the enumeration scheme of [22] for the nonrobust sparse portfolio problem: it can eliminate exponentially many subproblems in a single iteration and applies verbatim to that scheme. We introduce the sparsity-regularized robust portfolio problem in two alternative formulations and present the necessary background and notation in Section 2. Analytical results for both variants are presented in Sections 3 and 4, where we study structural properties of locally and globally optimal portfolios, including bounds on the number of nonzero entries and existence results. A branch-and-bound algorithm is proposed in Section 5. Finally, we compare our algorithm with a state-of-the-art solver for mixed-integer second-order cone programming problems in Section 6 and conclude in Section 7.
2
Problem Definition and Notation
Let IN = ({1, . . . , N }, <) be the strictly ordered index set, where < denotes the standard order. Any subset ω ⊆ IN inherits this property. We denote the all-ones vector by 1, the identity matrix by I, and the zero vector by 0. We denote the ith column of a matrix D by di . For ω ⊆ IN , the following notation for subvectors and submatrices will be used: rω := (r[ω[1]], . . . , r[ω[|ω|]]) ∈ R|ω| , Dω := (dω[1] )ω , . . . , (dω[|ω|] )ω ∈ R|ω|×|ω| . We define the zero-padding operator Zω : R|ω| → RN by ( 0, x = Zω (xω ), x[i] = xω [k],
i∈ / ω, if ω[k] = i.
We introduce the indicator ϕ : R → {0, 1} defined by ( X X 0, t = 0, ϕ(t) = so that ∥x∥0 := ϕ(x[i]) = ϕ(x[i]). 1, t ̸= 0, i∈I i∈σ(x) N
Here, ∥x∥0 = |σ(x)|, where σ(x) is the support of x (indices of nonzero entries), and |·| denotes P 1/p p cardinality. For x ∈ RN , the ℓp norm is defined for 1 ≤ p < ∞ as ∥x∥p := |x[i]| , i∈IN 3
∥x∥ √ ∞ = maxi∈IN {|x[i]|}, and the matrix induced norm by a positive definite matrix D is ∥x∥D = x⊤ Dx. Given ρ > 0, the open ℓp ball of radius ρ centered at x is Bp (x, ρ) := {y ∈ RN : ∥x − y∥p < ρ}. For a matrix A ∈ RM ×N , the spectral norm is ∥A∥2 = s1 (A), where si (A) is the ith singular value in decreasing order. We consider a financial market consisting of N risky assets and a single risk-free asset with deterministic period return rc . The vector of expected returns of the risky assets, denoted by r ∈ RN , is unknown, while their return covariance matrix D ∈ RN ×N is assumed to be known and positive definite. The initial wealth is normalized to one. Uncertainty in the mean returns is modeled through the ellipsoidal uncertainty set Ur̂ := r ∈ RN : ∥r − r̂∥D−1 ≤ γ , where r̂ denotes the nominal estimate of the mean return vector and γ > 0 is a tolerance parameter controlling the size of the uncertainty set. Let r̄ denote a prescribed target return level. The investor selects portfolio weights x ∈ RN for the risky assets and xc ∈ R for the risk-free asset so as to minimize portfolio variance while ensuring that the target return is achieved for all admissible realizations of the mean return vector. The resulting robust mean-variance portfolio optimization problem is formulated as min x⊤ Dx s.t. 1⊤ x + xc = 1 r⊤ x + rc xc ≥ r̄ ∀r ∈ Ur̂ (x, xc ) ∈ RN +1 . Under ellipsoidal uncertainty in the mean return vector, the semi-infinite robust constraint can be reduced to a single deterministic inequality. In particular, requiring that the portfolio achieves the target return for all admissible realizations of r, r⊤ x + rc (1 − 1⊤ x) ≥ r̄
∀r ∈ Ur̂ ,
is equivalent to enforcing that the nominal expected return, penalized by the worst-case deviation induced by the uncertainty set, exceeds the target level. This yields the deterministic robust counterpart r̂⊤ x + rc (1 − 1⊤ x) − γ ∥x∥D ≥ r̄, where the term γ ∥x∥D arises from minimizing the linear form r⊤ x over the ellipsoidal set Ur̂ . Consequently, the original semi-infinite robust mean-variance problem admits the equivalent finitedimensional formulation min x⊤ Dx
x∈RN
s.t. r̂⊤ x + rc (1 − 1⊤ x) − γ ∥x∥D ≥ r̄.
(1)
This equivalence follows from standard results in robust optimization, where linear constraints subject to ellipsoidal uncertainty admit exact deterministic robust counterparts involving norm penalties [6, 16]. Let β > 0 be a sparsity-inducing penalty parameter, and define the estimated excess return vector by r := r̂ − rc 1 together with the excess target return r̄ := r̄ − rc . We consider the following sparse robust mean-variance portfolio selection problem: min Fβ (x) := x⊤ Dx + β ∥x∥0
x∈RN
4
s.t.
r⊤ x − γ ∥x∥D ≥ r̄.
(P 1 )
An alternative robust formulation to the problem presented in (1) is obtained by reversing the roles of risk and return. Instead of minimizing variance subject to a robust return requirement, one may fix an admissible risk level and maximize the worst-case portfolio return. For a suitably chosen parameter T > 0, consider n o max min r⊤ x + (1 − 1⊤ x)rc s.t. x⊤ Dx ≤ T 2 . x∈RN r∈Ur̂
This formulation represents a robust counterpart of the classical mean-variance problem with a variance budget, and it provides an alternative scalarization of the risk-return trade-off. Exploiting the ellipsoidal structure of the uncertainty set Ur̂ , the inner minimization over r admits a closedform solution, yielding the deterministic equivalent max r̂⊤ x + (1 − 1⊤ x)rc − γ ∥x∥D
x∈RN
s.t.
x⊤ Dx ≤ T 2 .
Since r̂⊤ x + (1 − 1⊤ x)rc = r⊤ x + rc , the objective can be written in terms of excess returns as max r⊤ x + rc − γ ∥x∥D
x∈RN
s.t.
∥x∥D ≤ T,
and the constant rc can be dropped without affecting the set of optimal solutions. While this model captures the same robustness considerations as problem (1), the two are not equivalent. The original formulation enforces the target return as a hard robust constraint and minimizes risk accordingly, whereas the present model fixes the risk level and optimizes the worst-case return. Consequently, the two problems generally yield different optimal solutions, coinciding only for specific choices of T that recover the risk level induced by the optimal solution of problem (1). Using the same sparsity-inducing penalty parameter β as before, we introduce an alternative formulation aimed at promoting sparser solutions to the return maximization problem. Specifically, we consider the following optimization problem: min Gβ (x) := γ ∥x∥D − r⊤ x + β ∥x∥0
x∈RN
s.t.
∥x∥D ≤ T.
(P 2 )
In the following sections, we develop the theoretical foundations for both formulations. Our analysis covers local optimality conditions, rigorous bounds on the nonzero components of global minimizers, asymptotic results concerning the existence of global minimizers, and principled guidelines for selecting the parameter β to achieve prescribed sparsity levels.
3
Analysis of (P 1 )
We begin this section by introducing an assumption that rules out the degenerate solution x = 0, in which all wealth is invested in the risk-free asset. Assumption 1. The target return r̄ strictly exceeds the period return rc of the risk-free asset. We impose the following assumption to ensure that each asset’s nominal excess return is sufficiently large relative to its risk contribution and the level of mean uncertainty, thereby guaranteeing feasibility of the associated subproblems.
5
Assumption 2. The uncertainty radius γ satisfies |r[i]| > γ, min p i∈IN di [i] where
p
di [i] denotes the square root of ith diagonal entry of D.
Assumption 2 requires the magnitude of each asset’s nominal excess return relative to its volatility to exceed the uncertainty radius γ. Since short positions are permitted, this ensures that every singleton-support subproblem, and hence every nonempty support-restricted subproblem, is feasible. For ease of notation, define the feasible set of (P 1 ) as n o X = x ∈ RN : r⊤ x − γ ∥x∥D ≥ r̄ . (2) √ We also define H := r⊤ D−1 r, which represents the maximal achievable Sharpe ratio in the market and coincides with the slope of the capital market line. Remark 1. Note that for x ∈ RN and ω ⊇ σ(x), we have x⊤ Dx = x⊤ ω Dω xω . This property motivates searching for local minimizers by restricting to supports. For ω ⊆ IN , define Kω = {x ∈ RN : x[i] = 0, ∀i ∈ ω c }. We study the following subproblem to characterize local minimizers of (P 1 ): min x⊤ Dx.
x∈X ∩Kω
(Pω1 )
Since D is positive definite, the objective is strictly convex, and hence the problem has a unique solution whenever feasible.
3.1
Minimizers of (Pω1 )
We define the restricted feasibility set for a fixed ω ⊆ IN : n o Xω = u ∈ R|ω| : r⊤ ω u − γ ∥u∥Dω ≥ r̄ . Using these sets and the zero-padding operator, we define the equivalent problem to (Pω1 ) with a convex feasible region min u⊤ Dω u,
u∈Xω
|ω| ≥ 1.
(ZPω1 )
Remark 2. Assumption 2 is inherited by all subproblems:pfor any nonempty ω ⊆ IN , Dω is positive definite (as a principal submatrix of D) and |rω [k]| > γ Dω [k, k] holds for all k ∈ I|ω| , since the diagonal of Dω consists of diagonal entries of D. The following lemma and remark are crucial in establishing that Assumption 2 ensures the existence of feasible solutions for every subproblem associated with a subset ω ⊆ IN . Lemma 1. For any v ∈ RN , we have v ⊤ D−1 v = maxu∈RN {2u⊤ v − u⊤ Du}. Proof. Let v ∈ RN be fixed, and define Q(u) = 2u⊤ v − u⊤ Du. Then Q(u) = −(u − D−1 v)⊤ D(u − D−1 v) + v ⊤ D−1 v. Since D is positive definite, the first term is non-positive for every u ∈ RN , with equality if and only if u = D−1 v. Therefore the maximum is attained at u = D−1 v, and the desired identity follows. 6
Remark 3. Let u ∈ Xω , and write u = tv with t > 0 and ∥v∥Dω = 1. Then t(γ − r⊤ ω v) + r̄ ≤ 0. ⊤ Since r̄ > 0, feasibility requires γ − rω v < 0. This holds whenever q ⊤ −1 max rω v = r⊤ ω (Dω ) rω =: Hω > γ, ∥v∥Dω =1
p −1 where Hω = r⊤ ω (Dω ) rω is the Sharpe ratio of the subproblem defined by ω. To guarantee feasibility across all supports, it suffices to ensure minω̸=∅ Hω > γ. For any ω ̸= ∅ and i ∈ ω we use Lemma 1 by taking u = tei ∈ R|ω| and obtain: q p −1 Hω = r⊤ 2trω [i] − t2 di [i]. ω (Dω ) rω ≥ Maximizing the right hand side over t yields t = rω [i]/di [i] and we have q q |r [i]| |rω [i]| −1 r ≥ pω −1 r⊤ (D ) ⇒ r⊤ . ω ω ω ω (Dω ) rω ≥ max p i∈ω di [i] di [i] In fact, the minimum is attained on singletons: |r[i]| min Hω = Hmin := min p > γ. i∈I ω̸=∅ N di [i] Assumption 2 provides this bound. Thus, strict feasibility and convexity ensure a unique solution for (ZPω1 ). Proposition 1. For nonempty ω ⊆ IN , the unique solution of (ZPω1 ) is r̄ ξ(ω) := (Dω )−1 rω . Hω (Hω − γ) Proof. For the proof, see [23, Proposition 1]. Remark 4. For ω ⊆ IN with |ω| ≥ 1, we write Ξ(ω) = Zω (ξ(ω)) for its zero-padding to RN . The next proposition gives the optimal multiplier of the subproblem constraint and shows that the constraint is active at optimality. The multiplier is used in the branching rule of the algorithm. Proposition 2. For nonempty ω ⊆ IN , the optimal Lagrange multiplier associated with the constraint of (ZPω1 ) is µ∗ = (Hω2r̄−γ)2 . Moreover, the constraint is active at the optimal solution û, i.e., r⊤ ω û − γ ∥û∥Dω = r̄. Proof. The proof parallels the argument in [23, Proposition 1]. We also verify that the constraint 1 r⊤ ω u − γ ∥u∥Dω ≥ r̄ of (ZPω ) is active at optimality. Strict feasibility established above ensures that the KKT conditions apply. Suppose that the constraint is inactive at û. Complementary slackness then gives µ∗ = 0, and stationarity reduces to 2Dω û = 0. Since Dω ≻ 0, this implies û = 0, which is infeasible because r̄ > 0 by Assumption 1. Thus, the optimal solution û of (ZPω1 ) satisfies r⊤ ω û − γ ∥û∥Dω = r̄. Next, we show the dual multiplier in closed-form. By Proposition 1, we have Dω û =
r̄ rω , Hω (Hω − γ)
and
Dω û rω = . ∥û∥Dω Hω
Substituting these identities into the stationarity condition and rearranging the terms then yields µ∗ = (Hω2r̄−γ)2 . 7
The subproblem (ZPω1 ) is posed in the reduced space R|ω| . The following lemma confirms that solving it is equivalent to solving (Pω1 ) in RN . Lemma 2. Problems (ZPω1 ) and (Pω1 ) are equivalent. Proof. The map Zω : Xω → X ∩ Kω is a bijection. Moreover, for any xω ∈ Xω , we have x⊤ ω Dω xω = Zω (xω )⊤ DZω (xω ). Hence the two formulations are equivalent. Remark 5. For ω ⊆ IN with |ω| ≥ 1, the point Ξ(ω) ∈ RN is the unique solution of (Pω1 ).
3.2
(Local) Minimizers of (P 1 )
P Since Fβ (x) = x⊤ Dx + β i∈IN ϕ(x[i]), activating a zero component incurs an additional penalty β. The following proposition identifies a neighborhood in which this penalty cannot be offset by the change in the quadratic term of (P 1 ). Proposition 3. Let β > 0 and x̂ ∈ X . Define σ̂ = σ(x̂) and β ρ := min min |x̂[i]|, . i∈σ̂ 2(∥Dx̂∥1 + 1) Then ρ > 0, and: (i) If y ∈ B∞ (0, ρ), then
P
i∈IN ϕ(x̂[i] + y[i]) =
P
i∈σ̂ ϕ(x̂[i]) +
P
i∈σ̂ c ϕ(y[i]).
(ii) If y ∈ B∞ (0, ρ) ∩ (RN \ Kσ̂ ), then Fβ (x̂ + y) ≥ Fβ (x̂), with strict inequality whenever σ̂ c ̸= ∅. Proof. The proof follows identically from the argument in [22, Lemma 2]. Remark 6. Proposition 3 does not use feasibility of x̂ + y; in particular, it remains valid when y is restricted to perturbations with x̂ + y ∈ X . The following two results establish a correspondence between local minimizers of (P 1 ) and global minimizers of (Pω1 ) for ω ⊆ IN , in a manner analogous to Proposition 2 and Lemma 3 in [22]. The proofs are therefore omitted. Proposition 4. Let ω ⊆ IN , ω ̸= ∅. For any β > 0, the objective Fβ reaches a (local) minimum of (P 1 ) at Ξ(ω), with |σ(Ξ(ω))| ≥ 1 and σ(Ξ(ω)) ⊆ ω. Lemma 3. Let β > 0 and let x̂ be a (local) minimizer of (P 1 ). Then x̂ = Ξ(σ(x̂)). Proposition 4 and Lemma 3 together show that the set of local minimizers of (P 1 ) is precisely {Ξ(ω) : ∅ ̸= ω ⊆ IN }; in particular, every local minimizer is completely determined by its support.
3.3
Global Minimizers of (P 1 )
We begin by deriving a bounding box that contains all global minimizers of the problem. This box is subsequently employed for Big-M calibration and for tightening the feasible region. Proposition 5. Let β > 0 and suppose x̂ is a global minimizer of (P 1 ). Then we have r η r̄2 di [i] ∥x̂∥∞ ≤ < ∞, where η = min . 2 p i∈IN sN (D) |r[i]| − γ di [i] 8
Proof. Since x̂ is a global minimizer of (P 1 ), for any i ∈ IN , we have x̂⊤ Dx̂ + β ∥x̂∥0 ≤ Ξ({i})⊤ DΞ({i}) + β∥Ξ({i})∥0 ≤ Ξ({i})⊤ DΞ({i}) + β. Moreover, Assumption 1 implies r̄ > 0, hence 0 ∈ / X and therefore ∥x̂∥0 ≥ 1, which results in 2 2 ∥x̂∥D ≤ ∥Ξ({i})∥D for all i ∈ IN . By the definition of Ξ on singleton supports, Ξ({i}) =
r̄sign(r[i]) p ei , |r[i]| − γ di [i]
and therefore ( ∥x̂∥2D ≤ min ∥Ξ({i})∥2D = min i∈IN i∈IN
r̄2 di [i] p (|r[i]| − γ di [i])2
) = η.
For an upper bound on q the maximal entry, it is sufficient to rescale η with the smallest singular value, that is ∥x̂∥∞ ≤ sN η(D) . The primary objective of Proposition 5 was to identify a bounding box that contains all globally optimal solutions. Nevertheless, deriving an upper bound on the variance is of independent interest. The following result establishes a strict separation of nonzero components from zero under the prescribed parameters. The resulting lower bound on the nonzero entries of global minimizers constitutes the central insight underlying the proposed warm-start heuristic. Theorem 1. Let β > 0 and suppose x̂ is a global minimizer of (P 1 ). Let also σ̂ = σ(x̂) denote the support of x̂, and we define for all i ∈ IN ( ) ) ( p |r[i]| + γ di [i] r̄2 di [i] p p ρi = max ei + s ej and η = min . i∈IN s∈{−1,1} |r[j]| − γ dj [j] (|r[i]| − γ di [i])2 D j∈IN \{i}
Then for every i ∈ σ̂, we have ) (√ √ η+β− η r̄ p . , |x̂[i]| ≥ min ρi |r[i]| − γ di [i] Proof. Let i ∈ σ̂. We consider two cases based on the cardinality of σ̂. Case 1: Assume |σ̂| ≥ 2, then there exists j ∈ σ̂ such that i ̸= j. We define gij : RN → RN as gij (x) = x − x[i]ei − x[j]ej where ei and ej are canonical basis vectors of the specified index. We introduce the function f (ti , tj ) = Fβ (gij (x̂) + ti ei + tj ej ). Finally, we define a feasibility function h(ti , tj ) = r⊤ gij (x̂) + ti r[i] + tj r[j] − γ ∥gij (x̂) + ti ei + tj ej ∥D − r̄. Since x̂ is a global minimizer, Lemma 3 gives x̂ = Ξ(σ̂). By Remark 5, x̂ is therefore the unique solution of (Pω1 ) with ω = σ̂. Proposition 2 then implies that the restricted constraint is active at x̂. This, in turn, implies that the main constraint is also active, yielding h(x̂[i], x̂[j]) = 0. Precisely, we rewrite this equality as h(x̂[i], x̂[j]) = r⊤ x̂ − γ ∥x̂∥D − r̄ = 0 ⇒ r̄ = r⊤ x̂ − γ ∥x̂∥D .
9
Now, we define a map by fixing a choice for tj h(ti , tj ) = r⊤ gij (x̂) + ti r[i] + tj r[j] − γ ∥gij (x̂) + ti ei + tj ej ∥D − r̄ = γ ∥x̂∥D − ∥gij (x̂) + ti ei + tj ej ∥D + (ti − x̂[i])r[i] + (tj − x̂[j])r[j] ≥ −γ ∥(ti − x̂[i])ei + (tj − x̂[j])ej ∥D + (ti − x̂[i])r[i] + (tj − x̂[j])r[j] q p ≥ −γ |ti − x̂[i]| di [i] − γ |tj − x̂[j]| dj [j] + (ti − x̂[i])r[i] + (tj − x̂[j])r[j], where the first inequality exploits ∥u∥ − ∥v∥ ≥ −∥u − v∥ with u = x̂ and v p = gij (x̂) + ti ei + tj ej , and the second follows from the triangle inequality together with ∥ei ∥D = di [i]. Introduce the shorthand p |r[i]| + γ di [i] p δij := > 0, |r[j]| − γ dj [j] where positivity follows from Assumption 2. Set tj = sign(r[j]) |ti − x̂[i]| δij + x̂[j]. Substituting this choice into the lower bound above yields q p h(ti , tj ) ≥ −γ |ti − x̂[i]| di [i] + (ti − x̂[i])r[i] + δij |ti − x̂[i]| (|r[j]| − γ dj [j]) q p ≥ −(|r[i]| + γ di [i]) |ti − x̂[i]| + δij |ti − x̂[i]| (|r[j]| − γ dj [j]) = 0. Thus this choice ensures h(ti , tj ) ≥ 0. We may define a 1-dimensional restricted objective as f (t) = ∥gij (x̂) + tei + (sign(r[j]) |t − x̂[i]| δij + x̂[j])ej ∥2D + β ∥gij (x̂) + tei + (sign(r[j]) |t − x̂[i]| δij + x̂[j])ej ∥0 = ∥x̂ + (t − x̂[i])ei + sign(r[j]) |t − x̂[i]| δij ej ∥2D + β ∥x̂ + (t − x̂[i])ei + sign(r[j]) |t − x̂[i]| δij ej ∥0 . By the choice of tj above, the point x̂− x̂[i]ei +sign(r[j]) |x̂[i]| δij ej corresponding to t = 0 is feasible (h ≥ 0), while f (x̂[i]) = Fβ (x̂). Global optimality of x̂ therefore implies f (0) ≥ f (x̂[i]), hence we have ∥x̂∥2D + β ∥x̂∥0 ≤ ∥x̂ − x̂[i]ei + sign(r[j]) |x̂[i]| δij ej ∥2D + β ∥x̂ − x̂[i]ei + sign(r[j]) |x̂[i]| δij ej ∥0 . The vector x̂−x̂[i]ei +sign(r[j]) |x̂[i]| δij ej has at most |σ̂|−1 nonzero elements since the perturbation removes index i without adding a new nonzero index. Hence its ℓ0 -term is at most ∥x̂∥0 −1, yielding β ≤ ∥x̂ − x̂[i]ei + sign(r[j]) |x̂[i]| δij ej ∥2D − ∥x̂∥2D = (2x̂ − x̂[i]ei + sign(r[j]) |x̂[i]| δij ej )⊤ D(−x̂[i]ei + sign(r[j]) |x̂[i]| δij ej ) ≤ ∥2x̂ − x̂[i]ei + sign(r[j]) |x̂[i]| δij ej ∥D ∥−x̂[i]ei + sign(r[j]) |x̂[i]| δij ej ∥D ≤ 2 ∥x̂∥D + |x̂[i]| ∥ei − sign(x̂[i]r[j])δij ej ∥D |x̂[i]| ∥ei − sign(x̂[i]r[j])δij ej ∥D .
(3)
Here, the last two inequalities follow from the Cauchy-Schwarz and the triangle inequalities. By the √ definition of ρi , we have ∥ei −sign(x̂[i]r[j])δij ej ∥D ≤ ρi . Substituting this together with ∥x̂∥D ≤ η, established in the proof of Proposition 5, into (3), we obtain √ β ≤ (2 η + |x̂[i]| ρi ) |x̂[i]| ρi . 10
(4)
√ Consider the polynomial P (v) = ρ2i v 2 + 2 ηρi v − β. Since P is a quadratic with a positive leading coefficient, (4) implies that |x̂[i]| must be larger than the positive root of P . Evaluating these roots gives q √ √ √ −2 ηρi ± 4ηρ2i + 4ρ2i β − η± η+β = . v+ , v− = ρi 2ρ2i √
√ η+β− η
Hence, we obtain the lower bound |x̂[i]| ≥ . ρi Case 2: Assume |σ̂| = 1, then we have a closed form solution for the specific entry i ∈ σ̂, which takes the form x̂[i] =
r̄ sign(r[i]) p |r[i]| − γ di [i]
⇒
|x̂[i]| =
r̄ |r[i]| − γ
p . di [i]
In both cases the lower bound takes the form (√ ) √ η+β− η r̄ p , min , ρi |r[i]| − γ di [i] where the first term is active when |σ̂| ≥ 2 and the second when |σ̂| = 1. Since a global minimizer must fall into one of these two cases, the bound holds unconditionally for every i ∈ σ̂. √ The first term of the bound grows like β for large β, so the larger the sparsity penalty, the further the nonzero entries of a global minimizer with at least two assets must lie from zero. Having established bounds on the components of global minimizers, we next address their existence, for which we first verify that the objective is coercive. Lemma 4. The feasible set X defined in (2) is nonempty and closed, and its objective Fβ (x) is coercive on X . Proof. Since D ≻ 0, we have x⊤ Dx ≥ sN (D)∥x∥22 . Hence, Fβ (x) = x⊤ Dx + β∥x∥0 ≥ sN (D)∥x∥22 . Therefore, Fβ is coercive. Moreover, the feasible set X = {x ∈ RN : r⊤ x − γ∥x∥D ≥ r̄} is closed, since it is the preimage of a closed interval under a continuous function. Assumption 2 guarantees that X is nonempty. We are now ready to establish the existence of a global minimizer. Theorem 2. Let β > 0. Then the set of global minimizers of Fβ , X̂ = x̂ ∈ X : Fβ (x̂) = min Fβ (x) x∈X
is nonempty. Proof. By Lemma 4, the feasible set X is closed and nonempty, while Fβ is coercive. In addition, Fβ is lower semicontinuous, since ∥·∥0 is lower semicontinuous [21, Proof of Proposition 4.3]. Since Fβ is lower semicontinuous and coercive on the closed nonempty set X , it follows from [24, Theorem 1.9] that Fβ attains its minimum over X . Therefore, X̂ is nonempty. The following statement establishes the existence of a threshold value for the penalty parameter corresponding to each prescribed sparsity level. In particular, for every global minimizer of Fβ associated with a given sparsity, one can identify a penalty parameter that enforces that level. 11
Proposition 6. For any 1 ≤ k ≤ N − 1, there exists βk > 0 such that if β > βk , then every global minimizer x̂ of Fβ satisfies ∥x̂∥0 ≤ k. Proof. Fix k ∈ IN −1 and consider the set Xk+1 = {x ∈ X : ∥x∥0 ≥ k + 1}. Suppose first that Xk+1 is nonempty. Then, for each x ∈ Xk+1 , we have Fβ (x) = x⊤ Dx + β ∥x∥0 ≥ β(k + 1). Choose any support ω ⊂ IN such that |ω| ≤ k and Ξ(ω) ∈ X . Such a support always exists: by Assumption 2 and Remark 3, every singleton {i} with i ∈ IN satisfies H{i} = √|r[i]| > γ, which guarantees di [i]
Ξ({i}) ∈ X . Since |{i}| = 1 ≤ k for any k ≤ N − 1, taking ω = {i} for any i ∈ IN yields a valid choice. Choose βk so that βk ≥ Ξ(ω)⊤ DΞ(ω). For such a choice, Fβ (Ξ(ω)) = Ξ(ω)⊤ DΞ(ω) + β ∥Ξ(ω)∥0 ≤ βk + βk < β(k + 1) ≤ Fβ (x),
∀x ∈ Xk+1
whenever β > βk . Let x̂ ∈ X̂ denote a global minimizer of Fβ with x̂ ∈ X . Since Fβ (x̂) ≤ Fβ (Ξ(ω)) < Fβ (x),
∀x ∈ Xk+1 ,
we must have x̂ ∈ / Xk+1 . By the definition of Xk+1 , this implies ∥x̂∥0 ≤ k. If Xk+1 = ∅, then the existence of a global minimizer immediately yields ∥x̂∥0 ≤ k. Remark 7. For any β > 0, every admissible support induces a local minimizer of Fβ . Nevertheless, this correspondence is not one-to-one, as a single local minimizer may be generated by multiple supports via zero-padding operators. Consequently, the number of distinct local minimizers is bounded above by 2N − 1, which corresponds to the total number of nontrivial supports. In the next section, we extend the theoretical results developed here to the problem of robust return maximization under a variance budget. In particular, we adapt the main structural and sparsity-related properties to this risk-constrained setting and show that similar conclusions can be obtained in that framework.
4
Analysis of (P 2 )
In this section, we characterize the support-restricted, local, and global minimizers of (P 2 ) and derive bounds on the nonzero components of its globally optimal solutions. For ease of notation, define its feasible set as Y := x ∈ RN : ∥x∥D ≤ T , (5) We next impose a lower bound on T that rules out the zero vector as a global minimizer. Assumption 3. Under Assumption 2, we further assume ( ) p di [i]β p T > min . i∈IN |r[i]| − di [i]γ Remark 8. The purpose of this lower bound is to rule out the trivial solution x = 0, in which all wealth is invested in the risk-free asset. To see this, fix i ∈ IN and consider the singleton portfolio x=
T sign(r[i]) p ei . di [i]
12
This portfolio satisfies ∥x∥D = T , and its objective value is p T Gβ (x) = − p |r[i]| − γ di [i] + β < 0, di [i] where the strict inequality follows from Assumption 3. Thus, Gβ takes a negative value at a feasible point in Y, whereas Gβ (0) = 0. Hence, the zero vector cannot be globally optimal. Remark 9. The reasoning outlined in Remark 1 remains applicable to this objective function, as it preserves separability with respect to supports. This structural property plays a central role in developing a rigorous characterization of local minimizers. We consider the following restricted problem in order to characterize the local minimizers of (P 2 ): min γ ∥x∥D − r⊤ x.
x∈Y∩Kω
(Pω2 )
Under Assumption 2, we have Hω > γ for every nonempty ω. In this case, the subproblem admits a unique optimal solution, which we characterize in the next section.
4.1
Minimizers of (Pω2 )
For a fixed nonempty index set ω ⊆ IN , we introduce the corresponding restricted feasible region Yω = {u ∈ R|ω| : ∥u∥Dω ≤ T }. Based on this construction and the zero padding operator, we formulate a problem equivalent to (Pω2 ), now expressed over a convex feasible set: min γ ∥u∥Dω − r⊤ ω u,
u∈Yω
ω ̸= ∅.
(ZPω2 )
In this case, the subproblems (Pω2 ) and (ZPω2 ) admit a unique optimal solution, characterized below. Proposition 7. For nonempty ω ⊆ IN , the unique optimal solution of (ZPω2 ) is π(ω) :=
T (Dω )−1 rω . Hω
Proof. A proof of this result is given in [23, Proposition 2]. Remark 10. For ω ⊆ IN , we write Π(ω) = Zω (π(ω)) for its zero-padding to RN , similarly as in Remark 4. If Hω < γ, then u = 0 is the unique optimal solution of (ZPω2 ), that is, all wealth is held in the risk-free asset. Assumption 2 excludes this degenerate case, since it guarantees minω̸=∅ Hω > γ for every nonempty ω. As in the risk-minimization case (§3), the next proposition gives the optimal multiplier of the subproblem constraint and shows that the constraint is active at optimality. The multiplier is used in the branching rule of the algorithm. Proposition 8. For nonempty ω ⊆ IN the optimal Lagrange multiplier associated with the constraint of (ZPω2 ) is µ∗ = Hω − γ. Moreover, the constraint is active at the optimal solution û, i.e., ∥û∥Dω = T .
13
Proof. The argument follows closely that of [23, Proposition 2]. In addition, we establish that the restricted constraint ∥u∥Dω ≤ T of (ZPω2 ) is active at optimality. Suppose that the constraint is inactive at û. Complementary slackness then gives µ∗ = 0, so that û minimizes the unconstrained convex function u 7→ γ ∥u∥Dω − r⊤ ω u. If û ̸= 0, stationarity gives γDω û/ ∥û∥Dω = rω and taking −1 the Dω -norm of both sides yields Hω = ∥rω ∥(Dω )−1 = γ. If û = 0, the optimality condition rω ∈ γ ∂ ∥·∥Dω (0) gives Hω ≤ γ. Both conclusions contradict Hω > γ, which is guaranteed by Assumption 2. Consequently, the optimal solution û of (ZPω2 ) satisfies ∥û∥Dω = T . Dω û Stationarity at the nonzero optimal solution now gives rω = (γ+µ∗ ) ∥û∥ . Taking the Dω−1 -norm Dω
and using ∥û∥Dω = T yields Hω = γ + µ∗ . Hence, µ∗ = Hω − γ. Lemma 5. Problems (ZPω2 ) and (Pω2 ) are equivalent. Proof. The argument proceeds along the same lines as in Lemma 2.
4.2
(Local) Minimizers of (P 2 )
2 To analyze the local minimizers P of (P ), we express its objective using the indicator function ϕ; i.e., ⊤ Gβ (x) = γ ∥x∥D − r x + β i∈σ(x) ϕ(x[i]). The following proposition, analogous to Proposition 3, identifies a neighborhood in which activating components outside the support of a feasible point cannot decrease the objective of (P 2 ).
Proposition 9. Let β > 0 and x̂ ∈ Y \ {0}. Define σ̂ = σ(x̂) and ( ) β p ρ := min min |x̂[i]|, . i∈σ̂ ∥r∥1 + γ N ∥D∥2 + 1 Then ρ > 0, and: (i) If y ∈ B∞ (0, ρ), then
P
i∈IN ϕ(x̂[i] + y[i]) =
P
i∈σ̂ ϕ(x̂[i]) +
P
i∈σ̂ c ϕ(y[i]).
(ii) If y ∈ B∞ (0, ρ) ∩ (RN \ Kσ̂ ), then Gβ (x̂ + y) ≥ Gβ (x̂), with strict inequality whenever σ̂ c ̸= ∅. Proof. We first prove (i). For y ∈ B∞ (0, ρ) we have ∥y∥∞ < mini∈σ̂ |x̂[i]|. This implies ϕ(x̂[i] + y[i]) = ϕ(x̂[i]) for i ∈ σ̂. For i ∈ σ̂ c we have ϕ(x̂[i] + y[i]) = ϕ(y[i]), which gives the desired result. We now prove (ii). Let y ∈ B∞ (0, ρ) \ Kσ̂ . Then Gβ (x̂ + y) = γ ∥x̂ + y∥D − r⊤ (x̂ + y) + β ∥x̂ + y∥0 = Gβ (x̂) − r⊤ y + γ(∥x̂ + y∥D − ∥x̂∥D ) + β
X
ϕ(y[i])
i∈σ̂ c
≥ Gβ (x̂) − r⊤ y − γ ∥y∥D + β ∥yσ̂c ∥0 q ≥ Gβ (x̂) − ∥y∥∞ (∥r∥1 + γ N ∥D∥2 ) + β ∥yσ̂c ∥0 . For σ̂ c = ∅, inequality is trivial. If not, ∥yσ̂c ∥0 ≥ 1, and the given radius provides the inequality. The next result establishes a correspondence, analogous to Proposition 4 and Lemma 3, between local minimizers of (P 2 ) and global minimizers of (Pω2 ) for a given ω ⊆ IN , and the proofs are omitted accordingly. Proposition 10. Let ω ⊆ IN , ω ̸= ∅. For any β > 0, the objective Gβ reaches a (local) minimum of (P 2 ) at Π(ω), with |σ(Π(ω))| ≥ 1 and σ(Π(ω)) ⊆ ω. Moreover, any nonzero (local) minimizer x̂ of (P 2 ) satisfies x̂ = Π(σ(x̂)). 14
4.3
Global Minimizers of (P 2 )
We begin with an observation concerning the existence of an upper bound on the components of globally optimal portfolios. Remark 11. Let β > 0, and let x̂ be a global minimizer of (P 2 ). In Proposition 5 we developed a non-trivial upper bound on the nonzero entries of a global minimizer of (P 1 ). For (P 2 ), the feasible region Y defined in (5) is compact. Thus, the corresponding bound is immediate. Indeed, every global minimizer x̂ satisfies ∥x̂∥D = T and ∥x̂∥∞ ≤ √ T < ∞. sN (D)
We proceed by establishing lower bounds on the nonzero components of a globally optimal portfolio. Theorem 3. Let β > 0 and suppose x̂ is a global minimizer of (P 2 ). If σ̂ = σ(x̂) denotes the support of x̂, then for every i ∈ σ̂, we have ( ) β T p |x̂[i]| ≥ min ,p . |r[i]| + H di [i] di [i] Proof. Let |σ̂| ≥ 2 and i ∈ σ̂. We define gi : RN → RN as gi (x) = x−x[i]ei where ei is the canonical basis vector of the specified index. We introduce the function f (ti ) = Gβ (gi (x̂) + ti ei ). Finally, we define a feasibility function h(ti ) = T − ∥gi (x̂) + ti ei ∥D . Since x̂ is a global minimizer, it should be a minimizer of (Pσ̂2 ). As shown in Proposition 8, the restricted constraint is active at x̂. This implies that the main constraint is also active, yielding h(x̂[i]) = 0. Then, we have the following three cases: Case I: h(0) ≥ 0, in this case we have ∥gi (x̂)∥D ≤ T and due to the global optimality of x̂ we have γ ∥gi (x̂)∥D − r⊤ gi (x̂) + β ∥gi (x̂)∥0 ≥ γ ∥x̂∥D − r⊤ x̂ + β ∥x̂∥0 . β These together imply |x̂[i]| ≥ |r[i]| . β . di [i]+|r[i]| T f (0) < f (x̂[i]) and h(0) < 0. Since gi (x̂) ̸= 0 we can define û = ∥gi (x̂)∥ gi (x̂). û is a
Case II: h(0) < 0 and f (0) ≥ f (x̂[i]). Then, similarly we have |x̂[i]| ≥ √ γ
Case III: feasible point so we require Gβ (x̂) ≤ Gβ (û). Then we have
γT − r⊤ x̂ + β ∥x̂∥0 ≤ γT − r⊤ û + β(∥x̂∥0 − 1)
D
⇒
β ≤ r⊤ (x̂ − û).
Substituting the decomposition x̂ = gi (x̂) + x̂[i]ei and the definition of û, we obtain: T T ⊤ gi (x̂) = r[i]x̂[i] + 1 − r⊤ gi (x̂) β≤r gi (x̂) + x̂[i]ei − ∥gi (x̂)∥D ∥gi (x̂)∥D = r[i]x̂[i] + (∥gi (x̂)∥D − T )
r⊤ gi (x̂) . ∥gi (x̂)∥D
Since h(0) < 0, we have ∥gi (x̂)∥D > T , so the coefficient (∥gi (x̂)∥D −T ) is positive. By the definition ⊤
gi (x̂) of the dual norm, we have ∥gr i (x̂)∥ ≤ H. Applying this upper bound yields: D
β ≤ r[i]x̂[i] + (∥gi (x̂)∥D − T )H.
15
Next, using the triangle inequality ∥gip (x̂)∥D ≤ ∥x̂∥D + ∥x̂[i]ei ∥D and noting that ∥x̂∥D = T , we have ∥gi (x̂)∥D − T ≤ ∥x̂[i]ei ∥D = |x̂[i]| di [i]. Substituting this back into the inequality and using r[i]x̂[i] ≤ |r[i]| |x̂[i]|: β ≤ |r[i]| |x̂[i]| + |x̂[i]|
p p di [i]H = |x̂[i]| |r[i]| + di [i]H ⇒ |x̂[i]| ≥
β p . H di [i] + |r[i]|
We observe that this lower bound is suitable for all three cases above. Finally, if |σ̂| = 1 then we have closed form solutions for the specific entry i ∈ σ̂ x̂[i] =
T sign(r[i]) p di [i]
⇒
T |x̂[i]| = p . di [i]
Combining the bounds from the |σ̂| ≥ 2 and |σ̂| = 1 cases yields the desired result. Remark 12. The lower bound in Theorem 3 is simpler than the one in Theorem 1: it is linear in β, does not involve the pairwise quantities ρi , and depends on the data only through r[i], di [i], H and T . In our experiments it also separated the nonzero components of global minimizers from zero more clearly, which is the property exploited by the warm-start heuristic of Section 5.1. Because the feasible set is compact, the existence of a global minimizer can be verified more straightforwardly than in the earlier case. We first establish that the current objective function is lower semi-continuous for any selection of problem parameters. Lemma 6. For β > 0, the robust return maximization objective Gβ is lower semi-continuous. Proof. The function ∥·∥0 : RN → R is lower semi-continuous [21, Proof of Proposition 4.3]. Consequently, the objective function Gβ (x) = γ ∥x∥D − r⊤ x + β ∥x∥0 is also lower semi-continuous, as it is a sum of lower semi-continuous functions. Theorem 4. Let β > 0. Then the set of global minimizers of Gβ , Ŷ = x̂ ∈ Y : Gβ (x̂) = min Gβ (x) x∈Y
is nonempty. Proof. This result follows directly from the extended form of the Weierstrass’ extreme value theorem: a lower semicontinuous function attains its minimum over a compact feasible set. Finally, before concluding the theoretical developments, we observe that Proposition 6 and Remark 7 extend directly to the robust return maximization problem. Consequently, for any prescribed sparsity level k, there exists a regularization parameter βk > 0 that ensures the desired sparsity level in globally optimal portfolios. Furthermore, globally optimal portfolios can be identified among 2N − 1 locally optimal candidates, where this count arises from considering only nontrivial support sets. In the next section, we provide a brief description of the proposed algorithm and the associated heuristic procedures.
16
5
Enumeration Based BnB Algorithm
In this section, we develop an enumeration-based branch-and-bound algorithm for (P 1 ) and (P 2 ), following the scheme of [22] for the nonrobust counterpart and adapting it to the robust setting through the bounds established in Sections 3 and 4. We further refine the scheme with an additional pruning step in the bounding of right nodes, which allows an entire subtree to be replaced by a single leaf evaluation. We denote by P ⊆ IN the set of candidate assets that remain, whether or not the heuristic is applied, and decompose the problem into subproblems (Pω1 ) over support subsets ω ⊆ P . Each subproblem characterizes the local minimizers supported on ω and corresponds to a node in the enumeration tree. At each node we solve the lower-dimensional equivalent subproblem (ZPω1 ) and maintain a five-tuple (x, lb, ub, P, S), whose components are as follows. • x represents the solution obtained by solving (ZPω1 ) for a nonempty support subset ω. If ω = ∅, then the node is pruned. • lb represents a lower score, which is less than or equal to the upper bound at each node. For ease of reference, it will be referred to as the lower bound for the remainder of the paper. • ub represents an upper bound on the optimal value obtained from the corresponding node. • P denotes the set of candidate assets whose inclusion is still undecided, from which the branching asset is selected. • S denotes the set of assets already fixed in the support at this node. Each node produced by the algorithm is inserted into a priority queue, with its priority determined by the previously defined lower bounds. At the beginning of the next iteration, the node with the highest priority (i.e., the lowest lower bound) is removed from the priority queue, and the values of x, lb, ub, P , and S are updated according to this node. If multiple nodes share the same lowest lower bound, the node that was added to the queue first is selected. The branching process then proceeds from this chosen node, allowing the algorithm to systematically explore the search space. In the following subsections, we outline each step of the algorithm and present Algorithm 1, which summarizes the entire procedure. All subroutines can be applied to (P 2 ) with minor modifications, and in the following sections we present computational results for both problems.
5.1
Warm-Start Heuristic
Computation time increases with larger N , because the algorithm may encounter many suboptimal solutions, and the number of such solutions affects efficiency. To address this, a warm-start heuristic is introduced to reduce memory usage and problem size. The method first solves the problem with the full support set and determines a conservative elimination level using componentwise lower bounds. Then, for each asset i ∈ P , we compute the deficit b[i] − |uP [i]| between the componentwise lower bound b[i] of Theorem 1 and the weight uP [i] that asset i receives in the solution uP = ξ(P ) on the full candidate set, and remove the assets with the largest deficits from the candidate support set, thereby reducing the dimension of the problem early in the process. For instance, if the goal is a 10-sparse solution among 60 variables, about 30% of variables might be eliminated rather than removing all but 10, in order to avoid excluding potentially optimal components. This approach speeds up computation but introduces a trade-off between solution quality and runtime, so the elimination level must be chosen carefully. Guidance for this choice can come from 17
sparsity information provided by Proposition 6. If no warm-start is applied, that is, if no support is eliminated, the method reduces to the full branch-and-bound algorithm. The validity of the associated bounding and pruning scheme follows from the generic enumeration argument in [22, Appendix B]; consequently, when the algorithm terminates, the returned incumbent is globally optimal up to the prescribed stopping tolerance. The heuristic is most beneficial when the true solution is sparse, which aligns with portfolio optimization practice, since sparse portfolios reduce transaction costs.
5.2
Branching
When a node is taken from the queue for examination, the first step is to determine the most promising asset by analyzing the gradient of the Lagrangian function of (ZPω1 ). This branching strategy evaluates the quality of the current solution and aims to improve it. If the selected asset set at a node is empty, the procedure is modified by choosing the index that minimizes the varianceto-return ratio. For branching, we represent the Lagrangian function L by associating a Lagrange multiplier with the constraint defining Xω and compute its gradient (∇L) using Proposition 2. We then identify the most promising candidate among the set of possible assets. If the set of selected assets is nonempty, i.e., S ̸= ∅, we select j ∈ arg maxi∈P |∇Li | . In this case, rather than using the covariance matrix DP , we extract the submatrix DP,S , whose rows correspond to indices in P and columns correspond to indices in S, ensuring dimensional compatibility in the gradient computation. If S = ∅, the gradient rule is unavailable, since the submatrix DP,S is empty. In this case, we apply an alternative branching rule based on the variance-to-return ratio of assets in P , and select j ∈ arg mini∈P diag(DP )[i]/rP [i]. This rule favors assets with relatively low variance and high return. In both cases, an asset j is selected from the candidate set P , and the algorithm generates two distinct nodes: a left node in which the chosen asset is included in the support S L ← S ∪ {j}, and a right node in which the asset is excluded from the support P R ← P \ {j}, S R ← S ∪ P R . As a result, at every iteration the algorithm identifies the asset that appears most influential for the portfolio, which helps accelerate convergence relative to a standard BnB procedure. The primary objective is to obtain solutions with as few nonzero components as possible, that is, with small support sets.
5.3
Bounding
At the start of each iteration, the algorithm removes from the queue the tuple with the highest priority and updates the corresponding values (x, lb, ub, P, S). For left nodes, an upper bound is obtained by solving the subproblem with fixed support L S and evaluating an upper estimate using the objective value of (P 1 ) at the resulting solution, ubL = ξ(S L )⊤ DS L ξ(S L ) + β|S L |. The lower bound of a left node is computed by adding the sparsity penalty β to the parent node’s lower bound, lbL = lb + β. Depending on the triviality of the support, the node is then added to the queue with the updated bounds. For right nodes, the upper bound is inherited from the parent node because the branching step excludes an asset from the candidate set without generating a new feasible solution with fixed support. The lower bound is obtained from the relaxed subproblem of type (ZPω1 ). At initialization, a global lower bound is obtained by solving the problem with ω = IN , since this formulation omits the sparsity penalty and considers all indices. Following [22], after excluding the branching asset,
18
the right child is assigned the bounds lbR = ξ(S R )⊤ DS R ξ(S R ) + β |S| .
ubR = ub,
(6)
In the scheme of [22], the right child is then enqueued. Here, before enqueuing the right child, we apply the following additional pruning step. We first check whether lbR + β ≥ ub∗ − ϵ, where ub∗ denote the objective value of the incumbent, that is, the best feasible solution found so far, and ϵ is the prescribed optimality tolerance. If this condition holds, then the nonterminal descendants of the right subtree are pruned. The following proposition justifies this additional pruning step. Proposition 11. Consider (P 1 ) and let a node be represented by the current support set S and candidate set P . Let j ∈ P , P R = P \ {j}, and S R = S ∪ P R . Then lbR as defined in (6) satisfies the following property: if lbR + β ≥ ub∗ − ϵ, then every feasible portfolio x ∈ X such that S ⊊ σ(x) ⊆ S R satisfies Fβ (x) = x⊤ Dx + β∥x∥0 ≥ ub∗ − ϵ. Consequently, every descendant support S̃ with S ⊊ S̃ ⊆ S R may be pruned. Proof. Let x ∈ X satisfy S ⊊ σ(x) ⊆ S R , and define S̃ = σ(x). Since S̃ ⊆ S R , all nonzero components of x lie in S R . Therefore, the restriction of x to the indices in S R is feasible for (ZPω1 ) with ω = S R . By optimality of ξ(S R ), we obtain x⊤ Dx ≥ ξ(S R )⊤ DS R ξ(S R ). Since S ⊊ S̃, we also have ∥x∥0 = |S̃| ≥ |S| + 1. Hence, combining the two inequalities, we get Fβ (x) = x⊤ Dx + β∥x∥0 ≥ ξ(S R )⊤ DS R ξ(S R ) + β(|S| + 1) = lbR + β. Thus, if lbR + β ≥ ub∗ − ϵ, then Fβ (x) ≥ ub∗ − ϵ for every feasible x ∈ X with S ⊊ σ(x) ⊆ S R . Therefore, no feasible portfolio associated with a descendant support S̃ satisfying S ⊊ S̃ ⊆ S R can improve the incumbent by more than ϵ. Consequently, every such descendant support S̃ may be pruned. The effect of this pruning step can also be quantified. Indeed, consider the right child defined by branching on j, namely the node with candidate set P R = P \ {j} and current support S. Under the standard enumeration scheme [22], this node would be enqueued and explored. Since its support remains S and only indices from P R remain available for future branching, every support generated in that subtree satisfies S ⊆ S̃ ⊆ S R . When lbR + β ≥ ub∗ − ϵ, Proposition 11 shows that all supports with S ⊊ S̃ ⊆ S R may be pruned. The number of such supports is |P R |
X |P R | R = 2|P | − 1. i i=1
Therefore, the additional pruning step may remove an exponentially large portion of the right subtree in a single iteration. An analogous statement holds for (P 2 ) after replacing x⊤ Dx by γ∥x∥D − r⊤ x and ξ by π. The details are omitted for brevity.
19
Algorithm 1: A Branch-and-Bound Algorithm for (P 1 ) Input: tolerance ϵ = 1E − 8, local and global Subroutine: Branch&Bound(q, x, lb, ub, P, S) upper bounds ub = ub∗ = ∞, candidate support Input: queue q and node (x, lb, ub, P, S). set P = {1, . . . , N }, current support set S = ∅, if P ̸= ∅ then Branching: fixed-support solution uP = ξ(P ), lower bound 1 if |S| ≥ 1 then lb from solving (ZPω ) with ω = P , the corSelect j ∈ argmaxi∈P |∇Li |. responding componentwise lower bound b from else Theorem 1, and queue q = ∅. P )[i] Select j ∈ arg mini∈P diag(D . Warm-Start: let Jk ⊆ P be the set of inrP [i] dices corresponding to the k largest values of end b[i] − |uP [i]|, and update P ← P \ Jk . Enqueue Left node: S L ← S ∪ {j}. (uP , lb, ub, P, S) into q. Right node: P R ← P \{j}, S R ← P R ∪ S. ∗ while q ̸= ∅ and ub − lb > ϵ do Bounding the Right Node: Extract the highest-priority tuple if S R ≥ 1 then (x, lb, ub, P, S) from q. lbR = ξ(S R )⊤ DS R ξ(S R ) + β |S|. ∗ if ub < ub then if lbR + β < ub∗ − ϵ then Set ub∗ ← ub and x∗ ← x. Enqueue (x, lbR , ub, P R , S) into q. end end q ← Branch&Bound(q, x, lb, ub, P, S). end end Bounding the Left Node: return: ub∗ and its corresponding solution x∗ . if S L ≥ 1 then ubL = ξ(S L )⊤ DS L ξ(S L ) + β S L . Enqueue (ξ(S L ), lb + β, ubL, P R, S L ) into q. end end return: updated queue q.
20
5.4
Termination
The algorithm terminates when there are no remaining nodes to explore or when the best upper bound ub∗ found during the search matches the lower bound lb of the current node within a userdefined tolerance. In this case, further exploration cannot yield a better solution, and the procedure stops early.
6
Computational Results
The algorithms were evaluated using datasets sourced from [11] and [4]. The former includes daily price data (adjusted for dividends and stock splits) for the DowJones, EuroStoxx50, FTSE100, NASDAQ100, and S&P 500 indices. The latter comprises ETF, Eurobonds, and Italian Bonds datasets, providing daily asset returns derived from prices and total returns adjusted for dividends and splits. Summary statistics of these datasets are presented in Table 1. Experiments were run on a Linux cluster using a single CPU core and no GPU. The runs were carried out on nodes equipped with Intel Xeon Gold 6240 processors and 376 GiB RAM. The implementation used Python and Gurobi 13.0.2. The code used to generate all reported results is available at https://github.com/ecyayla/robust-mean-variance-portfolio-optimization. We benchmarked our method, enhanced with the warm-start heuristic, against Gurobi using the mixed-integer second-order cone reformulations of (P 1 ) and (P 2 ). For both approaches the stopping tolerance was set to 10−8 : the relative MIP gap tolerance for Gurobi, and the tolerance ϵ in the termination test ub∗ − lb ≤ ϵ for the branch-and-bound algorithm. To ensure a fair comparison, the number of threads used by Gurobi was restricted to one and its time limit was set to 12 hours. In the benchmarking experiments reported in Tables 2–3, the return of the risk-free asset rc is set to 0.0002, the target return r̄ to 5% above rc , and T = 1 for (P 2 ). Our support-wise analysis relies on Assumption 2, which guarantees feasibility of every nonempty support-restricted subproblem 2 of (P 1 ) and excludes degenerate support solutions p of (P ). We fix γ = 0.001 and retain, before running either method, the assets satisfying |r[i]| / di [i] ≤ γ. The same retained investment universe is used by BnB and Gurobi, and the asset counts in Table 1 report the resulting dimensions. Tables 2–3 summarize the computational performance of the proposed BnB algorithm and Gurobi. The columns labeled “BnB CPU Time (s)” and “Gurobi CPU Time (s)” report total CPU times in seconds. Entries marked with a dash (“-”) indicate that Gurobi did not certify optimality within 12 hours; its incumbent at termination is used to compute the reported sparsity and error. The “Solution Sparsity” column denotes the number of nonzero entries in the reported solution; when a single value is listed, both methods return solutions with the same sparsity, whereas when two values are listed, the first corresponds to BnB and the second to Gurobi. The column “Drop Rate” represents the proportion of assets removed during the warm-start phase, with 0 indicating that no elimination is performed. The column labeled “Error (%)” reports the relative objective error 100(vBnB − vG )/ |vG |, where vBnB and vG denote the objective values of the solutions returned by BnB and by Gurobi at termination, respectively. For (P 2 ), the optimal value is sometimes negative, which is why the denominator carries an absolute value. Since Gurobi is run under a fixed time limit, this value should be interpreted relative to its best available feasible solution when optimality is not certified. A negative value therefore indicates that our BnB method obtained a strictly better objective than the solution Gurobi returned at the time limit. When Gurobi terminates before the time limit, it certifies global optimality up to its gap and feasibility tolerances; in that case, no method can produce a better objective value for the same instance beyond these tolerances. Accordingly, our method is not intended to improve upon such solutions, but rather to deliver high-quality solutions within reasonable computation times, particularly for larger instances. 21
Finally, the column “Node Reduction (%)” reports the percentage reduction in the number of nodes explored by the branch-and-bound algorithm due to the additional pruning step of Proposition 11, relative to the standard enumeration scheme of [22]; a larger value indicates that a greater portion of the search tree is eliminated. Table 2 presents the corresponding results for (P 1 ) across datasets of varying sizes. On the small-scale instances such as ItalianBonds, ETF, DowJones, and EuroStoxx50, both methods return identical solutions with no elimination (drop rate 0); although the runtimes are small in absolute terms, BnB is already consistently faster than Gurobi. As the problem size increases (FTSE100 and NASDAQ100), the impact of the warm-start heuristic becomes more evident: by eliminating poorly performing assets it substantially reduces the problem size, and the returned solutions remain optimal (zero error) across all drop rates, with higher drop rates yielding the largest speedups and smaller drop rates increasing computation time. On FTSE100, Gurobi requires long computation times or fails to finish within the 12-hour limit, whereas BnB returns the same optimal solutions in seconds to a few minutes at higher drop rates, its runtime growing as the drop rate decreases. A similar pattern holds on NASDAQ100, where BnB outperforms Gurobi across all values of β and drop rates while matching its optimal objective. For the largest dataset, S&P500, Gurobi does not certify a solution within the time limit, while BnB returns solutions much faster but with a nonzero objective error (roughly 6% to 14%); here smaller drop rates reduce the error at the expense of additional computation time. Overall, the results suggest that the warm-start elimination strategy removes poorly performing assets and reduces the problem size substantially, while preserving optimality on all instances except the largest, where a controllable trade-off between computation time and solution quality remains. Finally, the additional pruning step substantially reduces the search effort on every instance, removing between 20% and 80% of the nodes explored by the enumeration scheme of [22]. Table 3 reports the results for (P 2 ), which follow a pattern comparable to that of (P 1 ). On the smaller datasets (ItalianBonds, ETF, DowJones, and EuroStoxx50) no assets are eliminated and the two methods return nearly identical portfolios, with BnB several times faster than Gurobi. For EuroStoxx50 with β = 5 × 10−4 , the reported −0.74% error is due to numerical rounding rather than a genuine improvement over Gurobi. On the larger instances (FTSE100 and NASDAQ100), Gurobi no longer certifies optimality within the 12-hour limit, and the warm-start heuristic becomes decisive: higher drop rates yield larger speedups at a small cost in accuracy, whereas lower drop rates reduce the error at the expense of additional computation time. On FTSE100 the error even becomes negative at the lower drop rates, where BnB improves upon the solution Gurobi returns at the time limit (for β = 10−3 , from 1.84% at drop rate 0.8 to −1.39% at 0.6), while on NASDAQ100 it stays below 1.3%. For the largest dataset, S&P500, BnB is again substantially faster, although the accuracy gap widens to between roughly 11% and 24% and narrows as the drop rate is reduced. Overall, the warm-start elimination substantially reduces the computation time of the larger instances and usually produces good-quality solutions, at the cost of some accuracy loss when the elimination is too aggressive. The node reduction is more modest here, ranging from 0% to roughly 42%, and it is largest on the instances whose optimal solutions are sparsest. This is not surprising because the pruning condition lbR + β ≥ ub∗ − ϵ is harder to satisfy when the portfolios have larger supports. Indeed, the reported solutions of (P 2 ) contain 5–271 nonzero components, compared with only 1–8 for (P 1 ).
6.1
Out-of-Sample Performance
We compare the robust (γ > 0) and purely sparse (γ = 0) models out of sample. We run a rolling-window backtest on the 100 size and book-to-market-sorted assets from the Kenneth French 22
Table 1: Dataset characteristics Index ItalianBonds ETF DowJones EuroStoxx50 FTSE100 NASDAQ100 S&P500
Num. of assets 11 24 28 46 82 70 420
Days 1564 1042 4276 4276 4276 4276 4276
Time interval 1/2013–12/2018 1/2015–12/2018 10/2006–2/2023 10/2006–2/2023 10/2006–2/2023 10/2006–2/2023 10/2006–2/2023
Table 2: Computational results comparison between BnB and Gurobi for (P 1 ). Solution Drop BnB CPU Gurobi CPU Error Node Sparsity Rate Time (s) Time (s) (%) Reduction (%) ItalianBonds 5E-5 1 0 0.01 0.05 0 35.35 1E-5 2 0 0.03 0.06 0 20.22 ETF 1E-4 2 0 0.06 0.27 0 57.43 5E-5 3 0 0.12 0.57 0 48.02 DowJones 5E-4 1 0 0.01 0.26 0 65.45 1E-4 2 0 0.02 0.42 0 71.30 EuroStoxx50 5E-4 1 0 0.02 0.67 0 79.78 1E-4 3 0 0.52 3.28 0 73.20 NASDAQ100 1E-4 2 0.5 1.20 56.24 0 70.30 1E-4 2 0.4 2.87 56.24 0 75.07 5E-5 3 0.5 5.44 264.82 0 61.92 5E-5 3 0.4 17.23 264.82 0 68.95 1E-5 8 0.5 132.76 35365.32 0 42.39 1E-5 8 0.4 1261.74 35365.32 0 44.01 FTSE100 5E-5 3 0.7 0.24 365.39 0 47.97 5E-5 3 0.6 1.26 365.39 0 54.00 5E-5 3 0.5 3.02 365.39 0 57.60 1E-5 8 0.7 4.40 0 38.91 1E-5 8 0.6 66.27 0 41.48 1E-5 8 0.5 386.43 0 48.11 S&P500 5E-5 3 0.9 37.64 - 14.03 67.85 5E-5 3-8 0.85 377.10 - 6.68 74.52 1E-5 8-3 0.9 6412.85 - 5.82 49.17 Dataset
β
23
Table 3: Computational results comparison between BnB and Gurobi for (P 2 ). Solution Drop BnB CPU Gurobi CPU Error Node Sparsity Rate Time (s) Time (s) (%) Reduction (%) ItalianBonds 5E-3 5 0 0.04 0.11 0 0.00 1E-3 6 0 0.03 0.20 0 0.00 ETF 1E-3 14 0 0.62 5.98 0 15.04 5E-4 17 0 0.28 2.81 0 11.20 DowJones 1E-3 9 0 0.32 7.46 0 41.61 5E-4 12 0 0.34 6.03 0 30.53 EuroStoxx50 5E-4 14-15 0 20.50 171.00 -0.74 30.42 1E-4 26 0 15.16 138.85 0 19.52 NASDAQ100 1E-3 17 0.5 538.27 - 0.48 29.25 5E-4 27-31 0.5 66.36 - 1.22 17.98 5E-4 30-31 0.4 1283.38 - 0.59 20.25 FTSE100 1E-3 11-17 0.8 0.07 - 1.84 3.74 1E-3 14-17 0.7 2.82 - -0.10 13.78 1E-3 16-17 0.6 154.19 - -1.39 19.32 5E-4 13-28 0.8 0.02 - 6.45 0.00 5E-4 21-28 0.7 1.06 - 2.44 4.30 5E-4 23-28 0.6 20.56 - 0.61 10.46 S&P500 1E-3 22-37 0.9 657.65 - 14.41 25.49 5E-4 31-73 0.9 169.11 - 22.62 12.67 1E-4 58-222 0.85 8.55 - 24.25 0.27 1E-4 78-222 0.8 564.04 - 16.67 0.00 5E-5 80-271 0.8 9.79 - 19.90 0.03 5E-5 120-271 0.7 5448.01 - 10.94 0.00 Dataset
β
Data Library1 (daily value-weighted returns, July 2004 to June 2026). At each rebalancing date, we estimate D and r from an in-sample window of 252 trading days and evaluate the resulting allocation out of sample over the next 21 trading days. Since the estimates are p recomputed at each date, Assumption 2 need not hold for the whole universe: assets with |r[i]| / di [i] ≤ γ are removed from the candidate set before solving, and rebalancing dates at which no asset survives this screening are skipped. The same screening and the same set of rebalancing dates are used for the robust and the sparse model, so that the two are compared on identical out-of-sample periods. This leaves 215 rebalances over the sample period. The risk-free return is set to rc = 5 × 10−5 and the target return to r̄ = 1.05 rc . Both models are solved by the branch-and-bound algorithm using the warm-start elimination of Section 5.1 at a drop rate of 0.3. We report the annualized out-of-sample Sharpe ratio, the annualized out-of-sample mean excess return, and the average number of selected assets |σ|. Across the tested configurations the robust model attains higher out-of-sample Sharpe ratios and mean excess returns than the sparse model in most settings. We do not claim that the robust model dominates the sparse one in general. Nonetheless, these results suggest that the robustness term can improve out-of-sample performance, which makes robust sparsity a worthwhile alternative to purely sparse portfolio selection.
7
Conclusion
In this paper, we considered a mean–variance portfolio selection problem that accounts for uncertainty in expected returns through an ellipsoidal uncertainty set, while at the same time promoting sparsity via an ℓ0 -penalty. Bringing these two aspects together leads to a problem that is both 1
Data available at https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html.
24
Table 4: Out-of-sample robust (γ > 0) versus sparse (γ = 0) performance on the Fama–French 100 Size×Book-to-Market portfolios. Mean excess returns are annualized and expressed in percent; |σ| is the average number of selected assets. Sharpe Mean return (%) |σ| β γ Sparse Robust Sparse Robust Sparse Robust 10−6 0.10 0.3160 0.4716 0.0102 0.0945 1.06 2.95 10−6 0.15 0.4505 0.5721 0.0111 0.2905 1.00 2.94 −7 5 × 10 0.10 0.4308 0.4312 0.0133 0.0862 1.32 3.75 5 × 10−7 0.15 0.4987 0.5704 0.0119 0.2896 1.06 3.38 nonconvex and discontinuous, and therefore inherently difficult to solve. To better understand this structure, we carried out a detailed analysis of both local and global minimizers for the two formulations studied, namely the robust risk minimization and robust return maximization models. In doing so, we clarified how local minimizers relate to support-restricted subproblems, established existence results for global solutions, and derived explicit lower and upper bounds on their components. These structural results guided the design of a tailored branch-and-bound algorithm. In particular, the bounds we obtained proved useful not only for pruning the search space, but also for constructing effective warm-start solutions. This combination plays an important role in keeping the computational effort manageable, even though the underlying problem is combinatorial due to the ℓ0 term. Our computational study on real financial data suggests that the proposed approach performs competitively against general-purpose mixed-integer second-order cone programming solvers, and in many cases achieves better performance in terms of running time without sacrificing solution quality. At the same time, the portfolios obtained reflect the intended balance: they are both robust to estimation errors and sparse enough to be practically implementable. Overall, the paper offers a unified perspective on combining robustness and exact sparsity in portfolio optimization under ellipsoidal uncertainty. There are several natural directions for future work, including the use of different uncertainty sets, the incorporation of additional practical constraints such as transaction costs or turnover limits, and the extension to alternative risk measures within the same framework.
Acknowledgments Buse Şen would like to acknowledge support as a part of NCCR Automation, a National Centre of Competence in Research, funded by the Swiss National Science Foundation (grant number 51NF40 225155).
References [1] On the computation of the efficient frontier in advanced sparse portfolio optimization. 4OR, 24:1–33, 2026. [2] Deniz Akkaya and Mustafa Çelebi Pınar. Minimizers of sparsity regularized Huber loss function. Journal of Optimization Theory and Applications, 187(1):205–233, 2020. [3] Deniz Akkaya and Mustafa Çelebi Pınar. Minimizers of sparsity regularized least absolute deviations. Journal of Global Optimization, 2025.
25
[4] Çağın Ararat, Francesco Cesarone, Mustafa Çelebi Pınar, and Jacopo Maria Ricci. MAD risk parity portfolios. Annals of Operations Research, 336:899–924, 2024. [5] Aharon Ben-Tal and Arkadi Nemirovski. Robust convex optimization. Mathematics of Operations Research, 23(4):769–805, 1998. [6] Aharon Ben-Tal, Laurent El Ghaoui, and Arkadi Nemirovski. Robust Optimization. Princeton University Press, 2009. [7] Dimitris Bertsimas and Ryan Cory-Wright. A scalable algorithm for sparse portfolio selection. INFORMS Journal on Computing, 34(3):1489–1511, 2022. [8] Dimitris Bertsimas and Melvyn Sim. The price of robustness. Operations Research, 52:35–53, 02 2004. [9] P. Bonami and M. A. Lejeune. An exact solution approach for portfolio optimization problems under stochastic and integer constraints. Operations Research, 57(3):650–670, 2009. [10] Joshua Brodie, Ingrid Daubechies, Christine De Mol, Domenico Giannone, and Ignace Loris. Sparse and stable Markowitz portfolios. Proceedings of the National Academy of Sciences, 106(30):12267–12272, 2009. [11] Francesco Cesarone, Rosella Giacometti, Manuel L. Martino, and Fabio Tardella. A returndiversification approach to portfolio selection. Computational Management Science, 22(2), 2025. [12] T.-J. Chang, Nigel Meade, John E. Beasley, and Yazid M. Sharaiha. Heuristics for cardinality constrained portfolio optimisation. Computers & Operations Research, 27(13):1271–1302, 2000. [13] Erick Delage and Yinyu Ye. Distributionally robust optimization under moment uncertainty with application to data-driven problems. Operations Research, 58(3):595–612, 2010. [14] Laurent El Ghaoui, Francois Oustry, and Hervé Lebret. Robust solutions to uncertain semidefinite programs. SIAM Journal on Optimization, 9(1):33–52, 1998. [15] Joel Goh and Melvyn Sim. Distributionally robust optimization and its tractable approximations. Operations Research, 58:902–917, 2010. [16] Donald Goldfarb and Garud Iyengar. Robust portfolio selection problems. Mathematics of Operations Research, 28:1–38, 2003. [17] Ken Kobayashi, Yuichi Takano, and Kazuhide Nakata. Cardinality-constrained distributionally robust portfolio optimization. European Journal of Operational Research, 309(3):1173–1182, 2023. [18] Miguel Sousa Lobo, Maryam Fazel, and Stephen Boyd. Portfolio optimization with linear and fixed transaction costs. Annals of Operations Research, 152(1):341–365, 2007. [19] Harry Markowitz. Portfolio selection. Journal of Finance, 7(1):77–91, 1952. [20] Richard O. Michaud. The Markowitz optimization enigma: Is ‘optimized’ optimal? Financial Analysts Journal, 45(1):31–42, 1989. [21] Mila Nikolova. Description of the minimizers of least squares regularized with ℓ0 -norm. Uniqueness of the global minimizer. SIAM Journal on Imaging Sciences, 6(2):904–937, 2013. [22] Buse Şen, Deniz Akkaya, and Mustafa Çelebi Pinar. Sparsity penalized mean–variance portfolio selection: Analysis and computation. Mathematical Programming, 211:281–318, 2025. [23] Mustafa Ç Pınar. On robust mean-variance portfolios. Optimization, 65(5):1039–1048, 2016. [24] R. T. Rockafellar and R. J.-B. Wets. Variational Analysis. Springer, 1998. [25] Wolfram Wiesemann, Daniel Kuhn, and Melvyn Sim. Distributionally robust convex optimization. Operations Research, 62(6):1358–1376, 2014. [26] Zhongming Wu, Kexin Sun, Zhili Ge, Zhihua Allen-Zhao, and Tieyong Zeng. Sparse portfolio optimization via ℓ1 over ℓ2 regularization. European Journal of Operational Research, 319(3):820–833, 2024. [27] Hongxin Zhao, Yilun Jiang, and Yizhou Yang. Robust and sparse portfolio: Optimization models and algorithms. Mathematics, 11(24), 2023.
26