Sparse High-Order Portfolios via Proximal DCA and SCA

In this paper, we aim at solving the cardinality constrained high-order portfolio optimization, i.e., mean-variance-skewness-kurtosis model with cardinality constraint (MVSKC). Optimization for the MVSKC model is of great difficulty in two parts. One is that the objective function is non-convex, the other is the combinational nature of the cardinality constraint, leading to non-convexity as well dis-continuity. Based on the observation that cardinality constraint has the difference-of-convex (DC) property, we transform the cardinality constraint into a penalty term and then propose three algorithms including the proximal difference of convex algorithm (pDCA), pDCA with extrapolation (pDCAe) and the successive convex approximation (SCA) to handle the resulting penalized MVSK (PMVSK) formulation. Moreover, theoretical convergence results of these algorithms are established respectively. Numerical experiments on the real datasets demonstrate the superiority of our proposed methods in obtaining high utility and sparse solutions as well as efficiency in terms of time usage.


INTRODUCTION
Portfolio management is a fundamental and challenging task for investors. One significant progress was made by Markowitz, who developed the mean-variance (MV) framework [1]. In the MV framework, the investors' purpose is to maximize their expected profit (i.e., mean return rate, or first moment) and minimize the corresponding risk (i.e., variance of portfolio, or second moment). Due to transaction costs, budget constraints, or even mental costs, investors may pick only a small number of assets out of all possible candidates, thus resulting in the sparse portfolios [2]. Some optimization methods have been developed for the MV framework with cardinality constraint [3][4][5].
The MV framework is based on the assumption that the returns follow a Gaussian distribution or investors have a quadratic utility [6]. However, in real financial markets, the Gaussian distribution assumption is seldom satisfied because asset returns are usually asymmetric and heavy-tailed [7,8].
To overcome this drawback, several studies took the third and fourth moments (i.e., skewness and kurtosis) into consideration, as they can better capture the asymmetry and heavy-tailed properties of a distribution [9][10][11]. They then proposed the mean-variance-skewness-kurtosis (MVSK) framework, which aims to maximize the mean and skewness while minimizing the variance and kurtosis. Assuming that investors have constant relative risk aversion (CRRA) [12] preferences, the MVSK model is the fourth-order Taylor expansion of the expected utility. Although the MVSK framework is more accurate than the MV framework, it requires solving a nonconvex optimization problem, which leads to computational challenges.
Recently, with the development in the optimization field, some algorithms were proposed to solve the MVSK-related problems. A difference-of-convex algorithm (DCA) was developed in [13] for the MVSK framework. Further improvement based on the difference-of-convex-sums-of-squares (DC-SOS) was proposed in [14]. Quite recently, to approximate the non-convex objective function more tightly, two algorithms based on the majorization-minimization (MM) and successive convex approximation (SCA) were proposed in [15]. By contrast, the cardinality constrained MVSK (MVSKC) framework has not received much attention yet. The work [16] considered a bi-objective optimization problem based on the trade-off between expected utility and cardinality. It is shown that there are gains in terms of outof-sample certainty equivalent and Sharpe ratio for certain cardinality levels. However, they applied a derivative-free solver based on direct multi-search [17], which may be cumbersome and sensitive to the initial point.
Our contributions. In this paper, based on the available structures of the MVSKC model, we consider integrating the cardinality constraint into the objective function via a penalty technique, which leads to the penalized MVSK (PMVSK) problem. Further, we present three approaches for solving the PMVSK problem, namely proximal difference-of-convex algorithm (pDCA), pDCA with extrapolation (pDCAe), and the successive convex approximation (SCA). Theoretically, we establish the convergence of pDCA and SCA to a stationary point of the problem at hand. Numerically, we demonstrate through experiments on a real dataset that our methods are able to obtain sparse solutions with lower objective function values. It is worth mentioning that SCA outperforms pDCA and pDCAe in terms of running time, which indicates that it is more favorable for large-scale problems.

PROBLEM FORMULATION
Suppose that the returns of N assets are given byr ∈ R N and the portfolio weights are given by w ∈ R N . The expected return of this portfolio, i.e., the first moment, is φ 1 (w) = E(w r) = w µ, where µ = E(r) is the mean return vector. We denote r = r − µ as the centered returns. The q-th central moment of the portfolio return is E[(w r) q ], from which we have the following: • The second moment, a.k.a. variance, of the portfolio return is is the co-skewness matrix. • The fourth moment, a.k.a. kurtosis, of the portfolio return is is the co-kurtosis matrix. Based on the above definitions, the MVSKC problem [16] is given as follows: min (1) where λ 1 , λ 2 , λ 3 , λ 4 > 0 are parameters to balance the four moments of the portfolio return, w 0 is defined as the number of non-zero elements of w, k < N is an integer controlling the number of assets to be selected, and α > 0 is to bound each element of w.
It is easy to verify that f is non-convex and twice continuously differentiable. The gradient and Hessian of each term in f are given in the following lemma.

Lemma 2. Under the constraints in
When it comes to the cardinality constraint, inspired by its DC property [18] w is the largest-k norm (i.e., the sum of the k largest elemens in absolute value), we transform (1) into the following PMVSK problem by taking the cardinality constraint as a penalty term: Here, ρ > 0 is the penalty coefficient.

ALGORITHM DESIGN
In this section, we propose three algorithms-pDCA, pD-CAe, and SCA-to solve the PMVSK problem (3).

pDCA
The main idea of pDCA is to successively construct a global upper bound of the objective function in (3) by linearizing the concave part. In the j-th iteration, we solve the subproblem where s j is a subgradient of w [k] at w j and can be computed efficiently through the following two steps: 1) sort the elements of |w j | in decreasing order, i.e., By introducing a new variable u ∈ R N , problem (4) can be cast as a convex quadratic programming (QP) problem: which can be efficiently solved by quadprog in MATLAB. In the rest of the paper, we will always utilize this technique to cast the 1 norm into linear inequality constraints. The complete description of pDCA is summarized in Algorithm 1.

pDCAe
Despite the common use of pDCA in many applications, it can be slow in practice [19]. To accelerate pDCA without increasing too much the computational cost, we adopt an extrapolation technique similar to that in FISTA [20] when approximating the concave part of the objective function. In the j-th iteration, we solve the subproblem where s j is a subgradient of w [k] at w j . The difference between (4) and (6) is that the latter linearly approximates the smooth concave part at y j , which is a carefully designed extrapolation point of w j−1 and w j . Note that a direct application of the extrapolation technique proposed in [19] for solving (3) can only lead to pDCA due to the nonconvexity of f . We summarize the pDCAe in Algorithm 2.

SCA
In this subsection, instead of utilizing the DC decomposition as in pDCA and pDCAe, we construct a strongly convex function that is not necessarily a global upper bound of f ncvx (w) but can lead to a tighter approximation. This is a Set y j = w j + β j (w j − w j−1 ).

3:
Perform the backtracking line search (8) to obtain the stepsize γ j .

7:
break 8: end if 9: end for kind of successive convex approximation strategy, and we refer to it as SCA for brevity.
Specifically, we havẽ and τ w > 0 guarantees thatf ncvx (w, w j ) is strongly convex. In addition to constructing a strongly convex local upper bound for f ncvx (w), we linearize the concave part −ρ w [k] . In the j-th iteration, we solve the following convex QP problem to getŵ j+1 : min Note that since problem (3) is non-convex and non-smooth, the traditional choice of stepsize [22] does not work. Hence, we perform a line search to guarantee the convergence of SCA. Specifically, given scalars 0 < c < 1 and 0 < β < 1, the stepsize γ j is set to be γ j = β mj , where m j is the smallest nonnegative integer m satisfying It is guaranteed that the stepsize γ j determined by (8) is nonzero and satisfies f p (w j+1 ) < f p (w j ) [23, Proposition 2]. The complete SCA algorithm is given in Algorithm 3.

THEORETICAL ANALYSIS
In this section, we provide convergence guarantees for our proposed pDCA and SCA in the following theorems. Theorem 1. [19, Theorem 4.1, Proposition 4.1] Let {w j } be the sequence generated by Algorithm 1. Then the following statements hold: • lim j→∞ w j+1 − w j = 0.
• Every limit point of {w j } is a stationary point of (3).
where Ω is the set of limit points of {w j }. Theorem 2. [23, Theorem 3] Let {w j } be the sequence generated by Algorithm 3. Then, every limit point of {w j } is a stationary point of (3).
We first show the convergence curves w.r.t. the iteration number and CPU time. Fig. 1 shows the objective value curves of pDCA, pDCAe, and SCA with ξ = 10 and the same random initialization. Fig. 2 presents the corresponding results when ξ = 5. Here, we only show the first 150 and 40 iterations of the algorithms for ξ = 10 and ξ = 5, respectively. We can observe from Fig. 1 (a) and Fig. 2 (a) that SCA converges fastest among the three algorithms, and pDCAe converges faster than pDCA. Since the QP subproblems (5)- (7), which constitute the main computational burden of the corresponding algorithms, require roughly the same cost, fewer iterations indicate less running time. This is verified in Fig. 1  (b) and Fig. 2 (b), which show that SCA takes the least time among the three methods.
Next, we compare our proposed algorithms with several baselines. One commonly used strategy, which we denote as relaxed MVSKC (RMVSKC), is to first relax the cardinality constraint to 1 norm constraint [15] (it is worth emphasizing that this kind of relaxation is not able to get a sparse solution) and then project the resulting solutionw * to satisfy the original constraints of MVSKC (1). The projection step amounts to solving To this end, we introduce a binary variable u to rewrite the cardinality constraint, leading to −αu ≤ w ≤ αu, and 1 u ≤ k. However, the resulting mixed-integer quadratic programming problem is challenging even with Gurobi (e.g., it takes more than 5 hours using Gurobi withw * ∈ R 50 ). Hence, we use the genetic algorithm (GA) instead. Utilizing the same transformation technique, we can also solve the MVSKC (1) via GA. From Table 1, we can see that pDCA, pDCAe, and SCA are superior in terms of the objective values obtained. In addition, it is achieved with much less time when using SCA. Table 1. f value and time usage of different methods with ξ = 10 and ξ = 5. Each method is tested five times using different initializations and the results are then averaged. 6. CONCLUSION In this paper, we considered the problem of high-order portfolio optimization with cardinality constraint. We proposed to recast the cardinality constraint into a penalized term and then developed three methods, namely pDCA, pDCAe, and SCA. Convergence results were established. Extensive experiments showed that our methods got lower objective values and are more efficient than the baselines.