domain: math math.math-oc info

text:
Massimiliano Kaucic
email: massimiliano.kaucic@deams.units.it

Filippo Piccotto
email: filippo.piccotto@phd.units.it

Gabriele Sbaiz
email: gabriele.sbaiz@phd.units.it

Giorgio Valentinuz
email: giorgio.valentinuz@deams.units.it

A hybrid level-based learning swarm algorithm with mutation operator for solving large-scale cardinality-constrained portfolio optimization problems

Keywords: level-based learning swarm optimizer, projection operator, mutation, exact penalty function, large-scale portfolio optimization

In this work, we propose a hybrid variant of the level-based learning swarm optimizer (LLSO) for solving large-scale portfolio optimization problems. Our goal is to maximize a modified formulation of the Sharpe ratio subject to cardinality, box and budget constraints. The algorithm involves a projection operator to deal with these three constraints simultaneously and we implicitly control transaction costs thanks to a rebalancing constraint. We also introduce a suitable exact penalty function to manage the turnover constraint. In addition, we develop an ad hoc mutation operator to modify candidate exemplars in the highest level of the swarm. The experimental results, using three large-scale data sets, show that the inclusion of this procedure improves the accuracy of the solutions. Then, a comparison with other variants of the LLSO algorithm and two state-of-the-art swarm optimizers points out the outstanding performance of the proposed solver in terms of exploration capabilities and solution quality. Finally, we assess the profitability of the portfolio allocation strategy in the last five years using an investible pool of 1119 constituents from the MSCI World Index.

Introduction

In modern portfolio theory, the classical mean-variance portfolio selection problem developed by Markowitz [START_REF] Markowitz | Portfolio selection[END_REF] plays a crucial role. Following this approach, investors should consider together return and risk, distributing the capital among alternative securities based on their return-risk trade-off. Since the pioneering work by Markowitz, the mean-variance optimization model (MVO) has been recognized as a practical tool to tackle portfolio optimization problems, and a large number of developments of the basic model has been investigated (see, for instance, [START_REF] Guerard | Handbook of portfolio construction. Contemporary applications of Markowitz techniques[END_REF] and [START_REF] Kolm | 60 years of portfolio optimization: practical challenges and current trends[END_REF]). Several authors have studied the multi-objective formulation of the MVO problem, in which the expected portfolio return is maximized and, at the same time, its variance is minimized ( [START_REF] Cesarone | A new method for mean-variance portfolio optimization with cardinality constraints[END_REF], [START_REF] Corazza | Particle swarm optimization with non-smooth penalty reformulation, for a complex portfolio selection problem[END_REF] and [START_REF] Kaucic | Equity portfolio management with cardinality constraints and risk parity control using multi-objective particle swarm optimization[END_REF]). These contributions aim to provide algorithms able to generate accurate dotted representations of all the sets of non-dominated portfolios in few iterations. On the contrary, in other studies, a single-objective formulation has been considered, maximising the return per unit of risk using a performance measure called the Sharpe ratio ( [START_REF] Sharpe | Mutual fund performance[END_REF], [START_REF] Farinelli | Beyond Sharpe ratio: optimal asset allocation using different performance ratios[END_REF] and [START_REF] Zakamouline | Portfolio performance evaluation with generalized Sharpe ratios: beyond the mean and variance[END_REF]). In this case, the attention is on the selection of the best-performing alternative in the set of feasible portfolios. However, from the point of view of a rational investor, the use of this indicator in periods of market downturns is questionable because it leads to prefer riskier portfolios ( [START_REF] Israelsen | A refinement to the Sharpe ratio and information ratio[END_REF] and [START_REF] Kaucic | Polynomial goal programming and particle swarm optimization for enhanced indexation[END_REF]).

In this paper, we focus on this occurrence and we propose a maximization problem where the objective function is a modified version of the Sharpe ratio, which is coherent with agent preferences also when risk premia are negative. We further consider four types of real-world constraints. First, a cardinality constraint manages the size of the portfolio; next, a budget constraint ensures that all the available capital is invested; lastly, bound constraints prescribe lower and upper bounds on the fraction of capital invested in each asset. To conclude, a turnover constraint implicitly controls the effect of the transaction costs on the portfolio rebalancing phases. We analyse the asset allocation problem from the perspective of an institutional investor who operates in large equity markets composed of hundreds or thousands of constituents and selects a restricted pool of stocks to build up a portfolio with a suitable performance with respect to the benchmark.

Since the resulting mixed-integer optimization problem has been proved to be NP-hard, finding possible optimal solutions becomes computationally challenging [START_REF] Moral-Escudero | Selection of optimal investment portfolios with cardinality constraints[END_REF]. For this reason, on the one hand, exact methods were proposed to supply optimal solutions, but they demand a significant amount of computation time when the problem size increases ( [START_REF] Li | Optimal lot solution to cardinality constrained mean-variance formulation for portfolio selection[END_REF], [START_REF] Shaw | Lagrangian relaxation procedure for cardinality-constrained portfolio optimization[END_REF] and [START_REF] Bertsimas | Algorithm for cardinality-constrained quadratic optimization[END_REF]). On the other hand, heuristic approaches can identify approximate and sometimes optimal solutions within reasonable computation time even when the problem size is huge ( [START_REF] Chang | Heuristics for cardinality constrained portfolio optimisation[END_REF], [START_REF] Crama | Simulated annealing for complex portfolio selection problems[END_REF], [START_REF] Woodside-Oriakhi | Heuristic algorithms for the cardinality constrained efficient frontier[END_REF] and [START_REF] Liagkouras | Examining the effect of different configuration issues of the multiobjective evolutionary algorithms on the efficient frontier formulation for the constrained portfolio optimization problem[END_REF]). In this context, swarm optimization algorithms, inspired by the self-organizing interaction among agents, have become popular in portfolio optimization theory in recent years [START_REF] Ertenlice | A survey of swarm intelligence for portfolio optimization: algorithms and applications[END_REF]. Specifically, the particle swarm optimization (PSO) algorithm has been widely employed to solve real-world financial problems since its first proposal ( [START_REF] Cura | Particle swarm optimization approach to portfolio optimization[END_REF], [START_REF] Zhu | A particle swarm optimization heuristic for the index tacking problem[END_REF], [START_REF] Thomaidis | A soft computing approach to enhanced indexation[END_REF], [START_REF] Zhu | Particle Swarm Optimization (PSO) for the constrained portfolio optimization problem[END_REF] and [START_REF] Corazza | A novel hybrid PSO-based metaheuristic for costly portfolio selection problems[END_REF]) due to its effectiveness in reaching optimal solutions. However, the aforementioned algorithm does not work efficiently when the problem size is large, leading to population stagnation and premature convergence [START_REF] Parsopoulos | On the computation of all global minimizers through particle swarm optimization[END_REF]. To improve PSO performance for large-scale optimization problems, several authors have designed many variants, such as competitive swarm optimizer [START_REF] Cheng | A competitive swarm optimizer for large scale optimization[END_REF], social learning particle swarm optimizer [START_REF] Cheng | A social learning particle swarm optimization algorithm for scalable optimization[END_REF] and level-based learning swarm optimizer (LLSO) [START_REF] Yang | A level-based learning swarm optimizer for large-scale optimization[END_REF]. In particular, the latter one has shown better exploitation ability in different environments. For this reason, we propose its use to solve our portfolio optimization problem. LLSO is inspired by the teaching concept that teachers should treat students differently according to their cognitive and learning abilities. Based on that, the general idea of the LLSO is to sort the swarm individuals in ascending order with respect to their fitness and then separate them into distinct levels. The best individuals are stored into the higher level and are not updated, preserving the most valuable information conveyed in the swarm. Unlike PSO, which uses the historically best positions to update the particles, LLSO employs predominant particles in the current swarm to guide the learning of the worst particles and to enhance the swarm diversity. Thus, particles in lower levels have more individuals in the upper levels to learn from and are focused on exploring the search space; those in higher levels mainly concentrate on the exploitation task. Even though LLSO shows promising capabilities in dealing with large-scale optimization problems, it is overly sensitive to its parameters. To mitigate this influence, an adaptive variant, henceforth ALLSO, has been introduced in [START_REF] Song | An adaptive level-based learning swarm optimizer for large-scale optimization[END_REF], which takes advantage of a swarm aggregation indicator to estimate the evolution state of the swarm. Two adaptive adjustment strategies are then applied to identify the best configuration setting for each generation.

Due to the fact that the swarm optimization algorithms are usually blind to the constraints, they have to be equipped with constraint-handling techniques [START_REF] Mezura-Montes | Constraint-handling in nature-inspired numerical optimization: past, present and future[END_REF] to be effective in real-world applications. A class of constraint-handling methods widely used in literature is represented by the penalty function methods, where a penalty term reduces the fitness value of the infeasible candidates. However, despite its simplicity, this method usually requires the definition of problem-dependent parameters that significantly impact algorithm performance. To overcome this issue, adaptive penalty techniques have been developed, in which the parameters are automatically set by using information gathered from the violated constraints at the current generation. For a more exhaustive overview of adaptive penalty techniques, we refer the reader to [START_REF] Barbosa | A critical review of adaptive penalty techniques in evolutionary computation[END_REF].

Therefore, to tackle the presented large-scale cardinality-constrained portfolio optimization problem, we combine the ALLSO with a novel hybrid constraint-handling technique, in which we integrate a projection operator into the self-adaptive penalty scheme developed by [START_REF] Costa | Theoretical and practical convergence of a self-adaptive penalty algorithm for constrained global optimization[END_REF]. Moreover, to further improve the exploitation power of the algorithm and the quality of the solutions, we introduce a novel mutation procedure, applied to the best individuals in the first level, which generalises the one inspected in [START_REF] Kaucic | Polynomial goal programming and particle swarm optimization for enhanced indexation[END_REF].

Even though similar asset allocation models with real-world constraints have been already studied in the literature, our work represents the first application of the LLSO algorithm for solving a largescale instance of the problem with the modified Sharpe ratio as the objective function. In addition, to the best of our knowledge, this is the first time such a hybrid constraint-handling technique and the mutation operator are presented and involved in the level-based learning paradigm.

Let us now give a more precise overview of the contents of the paper. The next section describes the investment framework, focusing on the portfolio optimization problem. In Section 3, we introduce the developed solver. More precisely, we first explain the adaptive LLSO and then we detail the proposed methods, namely the novel mutation operator and the hybrid constraint-handling technique. In the last part of the section, we summarise the entire procedure. In Section 4, we show the experimental results; in the last section, we depict the conclusions and future perspectives.

Portfolio design

Investment framework

Let us consider the standard portfolio selection problem introduced in [START_REF] Markowitz | Portfolio selection[END_REF]. We have a frictionless market in which no short selling is allowed, and all investors act as price takers. Assuming that n assets represent the investable universe, a portfolio is identified with the vector of assets weights x = (x 1 , . . . , x n ) ∈ R n , where x i ∈ R denotes the proportion of capital invested in asset i, with i = 1, . . . , n. Let R i be the random variable which stands for the rate of return of asset i, with expected value µ i . Hence, the random variable R p (x) = n i=1 R i x i indicates the rate of return of portfolio x. The expected rate of return of portfolio x is then defined as

µ p (x) = n i=1 x i µ i (2.1)
and its standard deviation, also called volatility, is given by

σ p (x) = n i=1 n j=1 c ij x i x j (2.2)
where (C) ij = c ij is the covariance between stocks i and j, with i, j = 1, . . . , n.

Since investors perceive large deviations from the portfolio mean value as damaging, (2.2) represents the so-called portfolio risk.

In such a setting, the portfolio choice is made only with respect to the expected portfolio rate of return and the portfolio risk, as stated in the following definition. Definition 2.1 Given two portfolios x, y, we say that x is preferred to y if and only if µ p (x) ≥ µ p (y) and σ p (x) ≤ σ p (y), with at least one strict inequality.

In other words, an investor prefers one portfolio to another if it has a higher expected rate of return and lower risk.

Objective function

In our portfolio selection problem, we take the Sharpe ratio ( [START_REF] Sharpe | Mutual fund performance[END_REF] and [START_REF] Sharpe | The Sharpe ratio[END_REF]) as an essential point of reference for measuring and comparing investment performance. It is defined as the ratio between the expected rate of return of the portfolio in excess of the risk-free rate and the standard deviation of the rates of return of the portfolio itself, that is

SR(x) = µ p (x) -r f σ p (x) (2.3)
where r f is the risk-free rate. This measure evaluates the compensation earned by the investor per unit of both systematic and idiosyncratic risks [START_REF] Caporin | A survey on the four families of performance measures[END_REF]. Thus, higher values of SR indicate more promising portfolios.

From a theoretical point of view, this choice is justified by the fact that several widely used performance measures are increasing functions of the Sharpe ratio ( [START_REF] Schuhmacher | Sufficient conditions for expected utility to imply drawdown-based performance rankings[END_REF] and [START_REF] Schuhmacher | A decision-theoretic foundation for reward-to-risk performance measures[END_REF]). Moreover, when the numerator in (2.3) is positive, this indicator is coherent with the risk-return profile of a rational investor. From a practical point of view, it can be easily calculated and its interpretation is simpler than most of recently proposed complex performance measures [START_REF] Auer | Performance hypothesis testing with the sharpe ratio: the case of hedge funds[END_REF].

However, as pointed out in [START_REF] Agudo | An analysis of Spanish investment fund performance: some considerations concerning Sharpe's ratio[END_REF] and [START_REF] Israelsen | A refinement to the Sharpe ratio and information ratio[END_REF], the reliability of this performance measure decreases when the excess rate of return is negative. In that case, one would prefer higher-risk portfolios using the Sharpe ratio. To overcome this issue, we adopt the following modification of (2.3), the so-called modified Sharpe ratio:

M SR(x) = µ p (x) -r f σ p (x) sign(µp(x)-r f ) (2.4)
where sign(z) is the sign function of z ∈ R. Observe that if the portfolio excess return is non-negative, the modified Sharpe ratio is equal to the Sharpe ratio. Otherwise, it multiplies the portfolio excess return by the standard deviation. In this manner, even in adverse conditions, portfolios with lower risk and higher excess return will be preferred.

Constraints

In our portfolio model, we consider the following constraints.

• Budget. All the available capital needs to be invested. In terms of portfolio weights, this translates to n i=1

x i = 1.

(2.5)

• Cardinality. We assume that the portfolio includes up to k assets, where k ≤ n. To model the inclusion or the exclusion of the i-th asset in the portfolio, a binary variable δ i is introduced as

δ i = 0, if asset i is excluded 1, if asset i is included (2.6)
for i = 1, . . . , n. The resulting vector of selected assets is δ = (δ 1 , . . . , δ n ) ∈ {0, 1} n , and the cardinality constraint can be written as

n i=1 δ i ≤ k.
(2.7)

• Box. A balanced portfolio should avoid extreme positions and foster diversification. Hence, we impose a maximum and a minimum limit for portfolio weights, that is

δ i l i ≤ x i ≤ δ i u i , i = 1, . . . , n (2.8)
where l i and u i are the lower and the upper bounds for the weight of the i-th asset, respectively, with 0 < l i < u i ≤ 1 to exclude short sales. • Turnover. To control the effect of the transaction costs in the portfolio rebalancing phases, we consider a portfolio turnover constraint. Let x 0 be a vector containing the current portfolio positions [START_REF] Shen | Doubly regularized portfolio with risk minimization[END_REF]. Then, the portfolio turnover constraint is

n i=1 |x i -x 0, i | ≤ T R (2.9)
where T R denotes the maximum turnover rate, which lies between 0 and 1. Note that if T R = 0 rebalancing is not allowed, and more trades are allowed when T R increases. where f (x) = -M SR(x).

3 Optimization algorithm

Adaptive level-based learning swarm optimizer

The algorithm evolves a swarm of N P candidate solutions using the so-called level-based population structure [START_REF] Yang | A level-based learning swarm optimizer for large-scale optimization[END_REF], according to which the evolution process is defined as follows.

1. At each iteration g, the individuals in the swarm are first sorted ascending based on their fitness and grouped into N L g levels, each one containing LP g = N P /NLg particles. In the last level, there are N P /NLg + N P %N L g particles. 1 Better individuals belong to higher levels, and a higher level corresponds to a smaller level index. Thus, L 1 represents the best level, while L N Lg is the worst one.

2.

To preserve the most valuable information conveyed in the current swarm, individuals belonging to L 1 are not updated and enter directly in the next generation. The p-th particle in level L l , denoted by x l, p (g), where l = 3, . . . , N L g and p = 1, . . . , LP g , is allowed to learn from two particles x l1, p1 (g), x l2, p2 (g) randomly extracted from two different higher levels L l1 and L l2 with l 1 < l 2 , and p 1 and p 2 are randomly chosen from {1, . . . , LP g }. For l = 2, we sample two particles from L 1 in such a way that x l1, p1 (g) is better than x l1, p2 (g) in terms of fitness function. Thus, the update rule for particle x l, p (g) is given component-wise by

v l, p i (g + 1) = r 1 v l, p i (g) + r 2 x l1, p1 i (g) -x l, p i (g) + φ g r 3 x l2, p2 i (g) -x l, p i (g) (3.1) x l, p i (g + 1) = x l, p i (g) + v l, p i (g) (3.2)
for i = 1, . . . , n, where v l, p i (g) denotes the i-th component of the velocity of particle p in level L l at generation g, and r 1 , r 2 , r 3 are real numbers randomly generated within [0, 1]. The parameter φ g ∈ [0, 1] controls the influence of the less performing exemplar x l2, p2 (g) on v l, p (g).

Based on [START_REF] Song | An adaptive level-based learning swarm optimizer for large-scale optimization[END_REF], both the parameters involved in the learning process at generation g, namely N L g and φ g , are adaptively adjusted based on the evolution state of the swarm by an aggregation indicator, which is defined as

s(g) = f g -f (x gbest (g)) f (x gbest (g)) + ξ (3.3)
where f g is the average fitness of the population at generation g, f (x gbest (g)) denotes the historically global best fitness up to iteration g, and ξ is a small positive value to avoid zero denominators.

Remark 3.1 When s(g) is high, particles are far from the current global best solution. Thus, the swarm is in an exploration phase. On the contrary, when s(g) is low, particles are close to the global best solution x gbest (g) and the swarm is in an exploitation phase.

To guarantee a control on the number of levels, N L g takes values in the set {N L min , . . . , N L max }, where N L min , N L max ∈ N are predefined lower and upper bounds. Moreover, to balance the level selection diversity and the exemplar diversity, N L g can be modified only when the relative improvement of the global fitness between generation g and generation g -1, given by

t(g) = f (x gbest (g -1)) -f (x gbest (g)) f (x gbest (g)) + ξ , ( 3.4)
slows down or stops, which corresponds to the cases t(g) < t(g -1) or t(g) = 0 respectively. The update of N L g then follows the rule

N L g = 2 • N L g-1 if s(g) < δ 1 2 • N L g-1 if s(g) ≥ δ (3.5)
where δ is a threshold in terms of the aggregation indicator to control the adjustment of N L g .

When N L g is out of the range, it is adjusted as follows

N L g =      N L rand if r < px N L max if r ≥ px and N L g > N L max N L min if r ≥ px and N L g < N L min (3.6)
where N L rand is uniformly sampled from {N L min , . . . , N L max }, r is a real number randomly generated within [0, 1], and px is a fixed probability employed to reset N L g . The update for φ g is designed in the following way

φ g = 0.35 + 0.1 • 1 1 + 10 • s(g) (3.7)
where s(g) is the value of aggregation indicator given in (3.3).

A preliminary numerical analysis reveals that a clamping procedure, limiting the magnitude of the velocity v l, p (g), provide a better exploration of the search space (in this regard, see also [START_REF] Oldewage | The merits of velocity clamping particle swarm optimisation in high dimensional spaces[END_REF]). This function can be written component-wise as

v l, p i (g) = min{max{v l, p i (g), v min i }, v max i }, (3.8)
where v min i and v max i are the minimum and the maximum velocity allowed for component i, with i = 1, . . . , n. In the experiments, recalling equation (2.8), we set

v max i = u i and v min i = -v max i .

Mutation operator

Instead of directly moving the individuals of the first level to the next generation, we propose to mutate them using an operator that combines two perturbation strategies properly developed for our portfolio optimization problem. More specifically, one technique is inspired by the swap operator proposed in [START_REF] Krink | Differential evolution and combinatorial search for constrained index-tracking[END_REF] and works as follows. First, we fix the maximum allowed number of non-null positions that could become zero, namely k swap max . Then, for each particle x 1, p (g) in level L 1 subject to swapping, we randomly sample from {1, . . . , k swap max } the number k swap of non-null positions that will be set to zero. At this point, for j = 1, . . . , k swap , let a j and b j be two randomly chosen positions in x 1, p (g), such that x 1, p aj (g) = 0 and l bj ≤ x 1, p bj (g) ≤ u bj . Thus, the modified individual, x1, p (g), is defined component-wise as

x1, p i (g) =          x 1, p i (g), if i = a j and i = b j l aj + x 1, p bj (g) -l bj u bj -l bj (u aj -l aj ), if i = a j 0, if i = b j .
(3.9)

In this paper, based on the preliminary experiments, k swap max = 0.05 • k , where k represents the maximum number of assets included in the portfolio.

Remark 3.2 This generalisation, allowing multiple swaps at the same time, improves the search capabilities of the original swap operator.

The other perturbation scheme focuses solely on the non-null components. For each x 1, p (g) to be mutated, let

I 1, p + (g) = {i : x 1, p i (g) > 0}
then, for all i ∈ I 1, p + (g), we define the interval

W 1, p i (g) = x 1, p i (g) -∆ i (g), x 1, p i (g) + ∆ i (g) (3.10)
where ∆(g) = 1 -g gmax+1 (ul), with g max be the maximum allowed number of iterations. The mutated component x1, p i (g) is randomly generated from the interval

W 1, p i (g) ∩ [l i , u i ]. For i / ∈ I 1, p + (g), we set x1, p i (g) = 0
. By narrowing the range over time, this procedure increases the exploration around the particles in L 1 .

Remark 3.3 By construction, the solutions modified by both the perturbation operators have at most k non-null positions and satisfy the box constraints.

For each particle in L 1 , the probability of applying the generalised swap operator decreases as the iteration counter increases according to the following rule

p swap (g) = 1 1 + exp(-0.005 • g) . (3.11)
In the initial stages, the proposed mutation favours the global search, using the generalised swap operator to identify the most promising subset of non-null decision variables. With the progress of the generations, the role of the refinement operator increases and, in the late stages, the algorithm focuses primarily on the local search.

The pseudo-code of the developed mutation procedure is reported in Algorithm 1.

Algorithm 1: Mutation procedure

Input : x 1, p (g), l, u, k swap max , g, ∆(g) Output: x1, p (g) 1 Set x1, p (g) = x 1, p (g) 2 Set I 0 = i : x 1, p i (g) = 0 3 Set I + = i : l i ≤ x 1, p i (g) ≤ u i 4 Calculate p swap (g) according to (3.11) 5 if rand() ≤ p swap (g) then 6 k swap → {1, . . . , k swap max } 7 for j = 1 to k swap do 8 a j → I 0 9 b j → I + 10 x1, p aj (g) = l aj + x 1, p bj (g) -l bj u bj -l bj (u aj -l aj )

Solution coding and hybrid constraint-handling procedure

Let us introduce some notation. Let C i denote a closed convex subset of R + , with i = 1, . . . , n, and K = {i 1 , . . . , i k } be any subset of indices of I = {1, . . . , n} with cardinality k ∈ N, so that I \ K is its complement in I. For all x ∈ R n , let x K be defined component-wise as

x K, i = x i , if i ∈ K 0, if i ∈ I \ K (3.12)
and let π K : R n → R k be the projection such that π K (x) = (x i1 , . . . , x i k ). We start by presenting the following proposition (in this regard, see also [START_REF] Zhang | A relaxed optimization approach for cardinality-constrained portfolios[END_REF]).

Proposition 3.4 Let y ∈ R n , with n ≥ 2.
Then, the optimal K for the problem

min x K : xi∈Ci 1 2 x K -y 2 (3.13)
is the set K * of indices corresponding to the k largest components of y.

The proof of this result is reported in Appendix A.

In other words, the proposition states that x K * is the vector with at most k non-null components which has minimum Euclidean distance from y among all x K , with K ⊂ I of cardinality k. Thanks to this projection, which implicitly enforces cardinality fulfillment, we can remove the vector of binary variables δ from the coding scheme of the solutions and we reformulate the portfolio optimization problem (2.11) only in terms of x K * . To this end, we introduce the set

B = x ∈ R n : x i = 0 or x i ∈ [l i , u i ] for i ∈ K * , x i = 0 for i ∈ I \ K * , n i=1 x i = 1 (3.14)
that is the set of the points satisfying all the constraints apart from the turnover condition. Further, let ψ(x) represent the value of the turnover function at x, which is given by

ψ(x) = n i=1 |x i -x 0,i | -T R. (3.15)
Then, the constrained optimization problem can be rewritten as

min x∈B f (x) s.t. ψ(x) ≤ 0. ( 3
δ * i = 1 if i ∈ K * 0 otherwise.
As previously observed, the standard ALLSO algorithm can only deal with unconstrained problems; thus, we propose to incorporate a hybrid constraint-handling technique in order to solve problem (3.16).

The building block of our procedure is based on the following lemma. The derivation of the orthogonal projection P [l, u] can be found in [START_REF] Beck | First-order methods in optimization[END_REF]. We now provide the main result concerning the projection phase. We refer the reader to Appendix A for the proof.

Proposition 3.7 Let y ∈ R n , with n ≥ 2, and K * * = {i ∈ K * : y i > 0}, with K * being the optimal set in Proposition 3.4. Assume that B in (3.14) is non-empty. Then, the orthogonal projection of y onto B is

P B (y) = π -1 K * * P [π K * * (l), π K * * (u)] (π K * * (y -η * 1) (3.18) where π -1 K (z) is the pre-image of z ∈ R |K * * | under π K and η * ∈ R is a solution of k i=1 P [π K * * (l), π K * * (u)] (y -η1) = 1. (3.19)
Let P g = {x p (g) ∈ R n : p = 1, . . . , N P } be the swarm at generation g, with g = 1, . . . , g max . Then, the proposed ALLSO variant maps the individuals in P g , which are updated using (3.1) and (3.2), onto the set B by means of the projector defined in (3.18). The resulting mutated swarm is denoted by Pg . Successively, we apply the self-adaptive penalty approach by [START_REF] Costa | Theoretical and practical convergence of a self-adaptive penalty algorithm for constrained global optimization[END_REF] to handle the turnover constraint and to guarantee the global optimality of solutions. More precisely, the objective function value at each projected individual in Pg , namely f (x p ), is normalized according to the formula

f (x p ) = f (x p ) -f min f max -f min where f min = min xp ∈ Pg f (x p ) and f max = max xp ∈ Pg f (x p ).
Similarly, the corresponding normalized constraint violation is given by

Ψ (x p ) =    max{ψ(x p ), 0} ψ max , if ψ max > 0 0, otherwise
where ψ max denotes the maximum of ψ(x p ) over all the mutated solutions in Pg which do not satisfy the turnover constraint. Finally, the penalty function is

F (x p ) =      f (x p ) if ψ(x p ) ≤ 0 f (ž) + R f Ψ (x p ) if ψ(x p ) > 0 and f (x p ) ≤ f (ž) f (x p ) + R f Ψ (x p ) if ψ(x p ) > 0 and f (x p ) > f (ž) , (3.20)
where R f represents the feasibility ratio for Pg , that is the percentage of individuals in Pg satisfying the turnover constraint. In (3.20), the reference point ž is a point belonging to Pg that satisfies the turnover constraint and has the lowest objective function value found so far. As in [START_REF] Costa | Theoretical and practical convergence of a self-adaptive penalty algorithm for constrained global optimization[END_REF], if the population has no feasible points, f (ž) is initially and temporarily set to f max , so that f (x p ) ≤ f (ž) for all xp ∈ Pg and f (ž) = 1. The value of f (ž) is updated only when the first feasible point is encountered.

We conclude this subsection by stating the following theorem, whose proof is omitted since it is similar to the one presented in [START_REF] Costa | Theoretical and practical convergence of a self-adaptive penalty algorithm for constrained global optimization[END_REF].

Proposition 3.8 The problem min

x∈B F (x)
with F as in (3.20), is equivalent to the problem (3.16).

Initialisation strategy and complete algorithm

For high dimensional problems, the common strategies of seeking a search space coverage by initializing the particles uniformly throughout the space as well as by increasing the size of the swarm are inefficient, because the search space grows exponentially with the dimension [START_REF] Van Zyl | A subspace-based method for PSO initialization[END_REF]. Moreover, the presence of highly constrained feasible regions in the search space exacerbates even more the initialization issue [START_REF] Datta | Evolutionary Constrained Optimization[END_REF].

To effectively address the low degree of feasibility in our portfolio rebalancing problem due to the complexity of the turnover constraint, we propose a direct initialization of the candidate solutions in a neighbourhood of x 0 . Let d min i and d max i be the minimum and the maximum allowed weight changes for x 0, i respectively, with i = 1, . . . , n. Let D p denote the total portfolio weight allowed to be re-allocated in x 0 for defining the p-th candidate solution x p (0), with p = 1, . . . , N P . Then, for each p, 1. we randomly select D p within [0, T R /2]; 2. we select a subset J -of k assets from the k assets with positive weight in x 0 , so that

x p j (0) = x 0, j -d j , for j ∈ J -
where d j is randomly sampled in d min j , d max j in such a way that j∈J -d j = D p , and x p j (0) = 0 or l j ≤ x p j (0) ≤ u j ; 3. we select a subset J + of k assets from the n -k assets with zero weight in x 0 , with k ≤ k , so that x p j (0) = x 0, j + d j , for j ∈ J + where d j is randomly sampled in d min , d max in such a way that j∈J + d j = D p , and l j ≤ x j ≤ u j ; 4. for j ∈ I \ (J -∪ J + ), we set x p j (0) = x 0, j .

The portfolios assembled using this scheme satisfy cardinality, box and turnover constraints. In this way, the initialization strategy encourages the swarm to focus on exploitation rather than exploration, thereby allowing it to identify promising solutions, even in problems with high dimension and small feasible regions.

Regarding the initial velocities, we set them all equal to the zero vector, that is v p (0) = 0, for p = 1, . . . , N P . The pseudocode of the proposed LLSO variant with adaptive parameters update, mutation of the particles in the first level and hybrid constraint-handling technique, shortly ALLSO-MUT-H, is reported in Algorithm 2. It can be noticed that, setting x 0 = 0 and T R = 1, ALLSO-MUT-H can also tackle portfolio optimization problems with no rebalancing. In this case, only the orthogonal projector is needed to move the unfeasible solutions to the feasible region.

Experimental analysis

This section is divided into two parts. We first point out the strengths and weaknesses of using the proposed algorithm to tackle large-scale cardinality-constrained portfolio optimization problems. The comparisons are made with two recent variants of the LLSO as well as with other state-of-the-art swarm optimization algorithms implementing the exact 1 -penalty function approach proposed in [START_REF] Corazza | A novel hybrid PSO-based metaheuristic for costly portfolio selection problems[END_REF]. Finally, we assess the profitability of the investment strategy in a real-world case study by varying the size of portfolios.

Algorithmic comparisons

For the algorithmic comparisons, we use three data sets from the OR-Library [START_REF] Canakgoz | Mixed-integer programming approaches for index tracking and enhanced indexation[END_REF], namely S&P 500, Russell 2000 and Russell 3000, which represent large capital market indices. In Table 1 are summarized the above cited data sets. For the calculation of the expected rates of return, we adopt a historical approach based on all the information available, that consists of 290 weekly prices for each asset. Then, to reduce the bias in the estimation of the covariance matrix C, we use the shrinkage estimator proposed in [START_REF] Ledoit | Honey, I shrunk the sample covariance matrix[END_REF].

Table 1: Data sets from [START_REF] Canakgoz | Mixed-integer programming approaches for index tracking and enhanced indexation[END_REF] with the corresponding number of weeks and number of market constituents (n) used in the estimation of parameters. For a fair of the solvers, all the algorithms have the same initial population of 500 individuals for each test set. In particular, we set d min i = 0.0005 and d max i = 0.0050 in the initialization strategy. For each algorithm, we perform 25 independent runs with 2000 generations. Moreover, all the portfolios from a given test set employ the following parameter setting. The risk-free value in (2.4) is set to zero. For the cardinality constraint (2.7), we consider k equal to 30% of the size n of the corresponding data set. The box thresholds in (2.8) are l i = 0.001 and u i = 0.05 for each asset. Regarding the turnover constraint (2.9), T R is set equal to 0.20 and the vector of current positions x 0 is fixed for all the compared algorithms and in all simulations by randomly sampling once from each set of feasible portfolios.

Mutation effects on LLSO variants

Our first task is to study the impact of the developed mutation operator on the LLSO-type algorithms, all equipped with our hybrid constraint-handling technique. For this purpose, we compare the following variants of the LLSO: the dynamic LLSO (DLLSO) [START_REF] Yang | A level-based learning swarm optimizer for large-scale optimization[END_REF], the adaptive LLSO (ALLSO) [START_REF] Song | An adaptive level-based learning swarm optimizer for large-scale optimization[END_REF] and the reinforcement learning level-based particle swarm optimization (RLLPSO) algorithm [START_REF] Wang | Sun: A reinforcement learning level-based particle swarm optimization algorithm for large-scale optimization[END_REF]. The parameters setting of each algorithm is chosen following the literature, as reported in Table 2.

Table 2: Parameter settings of the algorithms used in the comparisons. Table 3 shows four performance metrics linked to the best objective function value over the 25 runs, on the three public data sets. The best results are highlighted in bold font. We note that the ALLSO-MUT-H outperforms the competitors in all the case studies, presenting the lowest mean objective function value. Further, we remark that all the inspected LLSO variants are able to find feasible solutions. Focusing on the mutation benefits, we observe from Table 4 that the mutation has a significant impact on the performance of solvers. Specifically, although mutation-based optimizers present higher volatility than their counterparts, they always show lower results in terms of minimummaximum range of the best solutions.

)
Figure 1 shows the convergence and the diversity analyses of the compared solvers on the three data sets. From the first set of graphs, we note that the three mutated algorithms are able to reach significantly lower objective function values, and the ALLSO-MUT-H performs better than the others. Moreover, the algorithms without mutation show population stagnation around 100 generations, meaning that they converge to a local minimum and are not able to further explore the search space. This is confirmed by the results showed in the logarithmic scale plots of the diversity measures. We can observe that the ALLSO-MUT-H and the DLLSO-MUT-H are able to escape from the local minima, due to the oscillatory behaviour of the swarm diversity.

Comparison with state-of-the-art swarm optimization algorithms

In the previous subsection, we have analysed the impact of the mutation on the capabilities of LLSObased algorithms, finding that the ALLSO-MUT-H is the more efficient choice in terms of convergence and quality of solutions. Now, the aim is to compare the above quoted optimizer with other stateof-the-art swarm optimization algorithms, namely the PSO [START_REF] Corazza | A novel hybrid PSO-based metaheuristic for costly portfolio selection problems[END_REF] and the Firefly algorithm (FA) [START_REF] Yang | Firefly algorithm, stochastic test functions and design optimisation[END_REF], both endowed with an exact 1 -penalty function. In this regard, the literature presents a wide range of penalty methods to tackle constraint-handling problems [START_REF] Mezura-Montes | Constraint-handling in nature-inspired numerical optimization: past, present and future[END_REF]. However, in the context of portfolio optimization allocation, the technique proposed in [START_REF] Corazza | A novel hybrid PSO-based metaheuristic for costly portfolio selection problems[END_REF] is the only one for which the convergence is guaranteed. For this reason, we have decided to adopt the exact 1 -penalty technique.

We recall also that we adopt the same parameter setup presented above in Table 2. We exhibit the statistics of the comparison in Table 5, where are displayed the percentage of feasible solutions provided by the different solvers over the 25 runs; the mean of the constraint violation function CV for the non-feasible solutions; the average value of the penalty function F 1 over the 25 runs. Notice that, the penalty function corresponds to the objective function when the solutions are feasible. Looking at the results, we can argue that the ALLSO-MUT-H reaches the best mean value of the penalty function in all the data sets, and it always provides feasible solutions. This insight is confirmed by the convergence analysis plots in Figure 2, which show the benefits of our constraint-handling technique in terms of accuracy of solutions. Moreover, the diversity graphs suggest that our solver is the sole algorithm able to exhibit exploration and exploitation phases alternatively.

Real-world application

Data description and investment setting

The constituents of the MSCI World index, on 31st January 2022, form our investible pool. The data set has been downloaded from DataStream and consists of monthly prices covering the period from January 2012 to January 2022 for a total of 121 months. Stocks with missing observations were disregarded, and thus the final data set includes 1119 stocks. For the performance comparisons, we introduce a value-weighted benchmark index with the same constituents and, as in the previous analysis, we set r f = 0. The portfolio design employs the following parameter setting. For the cardinality constraint (2.7), we consider k ∈ {335, 167, 55, 22}, corresponding to k % = 30%, 15%, 5% and 2% of the pool size, respectively. As stated in the introductory part of Section 4.1, we recall that the box thresholds in (2.8) are l i = 0.001 and u i = 0.05 for each asset i, with i = 1, . . . , 1119, and the turnover rate in (2.9) is set equal to 0.20, as in [START_REF] Kaucic | Polynomial goal programming and particle swarm optimization for enhanced indexation[END_REF].

We use a rolling time window procedure to rebalance optimal portfolios every month, from January 2017 to January 2022, to point out the effects of the market changes on the behaviour of the investments and, as a consequence, the total number of ex-post dates is 61. We solve the corresponding problem instances by employing overlapping 60-months windows, which are updated every month by removing the oldest data and including the latest information.

In each quoted window, as already pointed out above, we adopt a historical approach to calculate expected rates of return, and to reduce the bias in the estimation of the covariance matrix C, we take advantage of the shrinkage estimator proposed in [START_REF] Ledoit | Honey, I shrunk the sample covariance matrix[END_REF]. Let us denote by x t the optimal portfolio at the ex-post month t, with t = 1, . . . , 61. Due to the time dependence of the considered investment plan, we rewrite the turnover constraint (2.9) as follows

n i=1 |x t,i -x t -,i | ≤ T R . (4.1)
In the previous equation, x t -= (x t -,1 , . . . , x t -,n ) represents the portfolio to be rebalanced [START_REF] Shen | Doubly regularized portfolio with risk minimization[END_REF], which is defined for t = 2, . . . , 61 as

x t -, i = x t-1, i R g t-1, i n j=1 x t-1, j R g t-1, j (4.2)
with the denominator being the gross portfolio return at month t -1.2 At time t = 1, we set x t -= 0 and T R = 1.

Let us assume a self-financing strategy with an initial wealth W 0 = 10, 000, 000 $. Then, we explicitly evaluate the magnitude of the trading through the cost function λ(x t , x t -) introduced in [START_REF] Beraldi | Dealing with complex transaction costs in portfolio management[END_REF]. As reported in Table 6, we consider the transaction cost structure characterized by decreasing cost rates as the traded value increases. Trading segment ($) Fixed fee ($) Proportional cost (%) 0 -7,999 40 0 8,000 -49,999 0 0.5 50,000 -99,999 0 0.4 100,000 -199,999 0 0.25 ≥ 200, 000 400 0

Ex-post performance measures

The following measures are considered to evaluate the profitability of the investment strategies. Let r out p,t be the ex-post portfolio rate of return realized at time t, with t = 1, . . . , 61. First, we consider the so-called ex-post Sharpe ratio [START_REF] Sharpe | The Sharpe ratio[END_REF], defined as

SR out = µ out σ out (4.3)
where µ out and σ out are the mean and the standard deviation of the ex-post portfolio rates of return, respectively.

The second measure employed in the analysis is the so-called Omega ratio [START_REF] Keating | A universal performance measure[END_REF], defined as the ratio between the gains over a threshold level and the losses under a threshold level. In this study, we set both thresholds equal to zero, that is The information gathered from these performance measures draws a complete picture of the ex-post portfolio return distribution. In particular, the ex-post Sharpe ratio describes the central part of the portfolio return distribution, while the Omega ratio considers the behaviour of profits and losses. Further, to measure the profitability of the investment at time t, we compute the net wealth as

Omega = 61 t=1 r out p,
W t = W t-1 1 + r out p,t -λ(x t , x t -). (4.5)
Then, we compare the profitability of the investments using the so-called compound annual growth rate [START_REF] Sharpe | The Sharpe ratio[END_REF], which in our case is calculated as

CAGR = W 61 W 0 12 61 -1 (4.6)
where W 0 represents the initial wealth and W 61 is the final wealth.

To evaluate the capacity of a strategy to avoid high losses, we introduce the drawdown measure [START_REF] Chekhlov | Drawdown measure in portfolio optimization[END_REF], which can be written

DD t = min 0, W t -W peak W peak (4.7)
where W peak is the maximum amount of wealth reached by the strategy until time t. In particular, we consider the mean and the standard deviation of the drawdown measure over time. Finally, we propose to measure the effect of the costs on the available capital in the out-of-sample period by

Λ % = 1 61 61 t=1 λ(x t , x t -) W t-1 • 100. (4.8)

Ex-post performance analysis

In the ex-post analysis, we investigate how the performance of the proposed asset allocation model changes by varying the k parameter in (2.7). First, we remark that, for any ex-post dates, the proposed hybrid LLSO variant identifies feasible solutions for all the portfolio sizes. The empirical results are summarised in Table 7, where the number of assets of the considered strategies is also displayed. Note that for each value of k, the proposed investments provide better performances than the value-weighted benchmark. This result implies that introducing a cardinality constraint in the portfolio model allows to choose a subset of the most profitable assets in the investible pool.

In terms of the return-risk profile, strategies with k % = 30%, 15%, and 5% show comparable performances, while the strategy with k % = 2% has a lower Sharpe ratio, which is due to its large volatility. Similar conclusions can be made about the Omega ratio, which expresses the gain-loss profile of the strategies. Despite better performance with respect to Sharpe and Omega ratios, strategies involving portfolios with a larger number of assets generate less wealth. Moreover, we observe that reducing cardinality leads to more profitable portfolio strategies.

Concerning the drawdown measures, the 5% asset allocation model is the most conservative, while the one with k % = 2% is the worst. Thus, the performance deteriorates by reducing portfolio size below a critical threshold.

As highlighted in the last two rows of Table 7 and in Figure 3, the impact of transaction costs for strategies with small k is negligible. On the contrary, portfolios with many assets have more fluctuations in the rebalancing phases, leading to higher trading commissions with a significant impact on the wealth generated.

Summing up, we can infer that the strategy with k % = 5% shows the best balance between riskadjusted performance measures and capability to generate net profits.

Conclusions and future works

In this paper, we have developed a swarm optimization algorithm for solving a large-scale cardinalityconstrained portfolio optimization problem, where the modified Sharpe ratio performance measure represents the objective function. We have considered four real-world constraints: cardinality, box, budget and turnover constraints. Due to the properties of the model inspected, we have proposed a variant of the LLSO equipped with a hybrid procedure to manage the constraints efficiently. Moreover, a novel mutation operator has been introduced to improve the accuracy of solutions. Our solver capabilities have been compared with those of two variants of the LLSO as well as other state-ofthe-art swarm optimization algorithms endowed with an 1 -penalty function. Numerical experiments on three publicly available large-scale data sets showed the outperformance of our hybrid procedure.

From the financial point of view, we have analysed the sensitivity of the portfolio model to the cardinality constraint with data from the last five years of the MSCI World index. We have found that portfolios of small size are more competitive with respect to the value-weighted benchmark index, also in periods of market downturns. Specifically, the losses are reduced, and the cost impact on the available capital is marginal compared to the profits. In the future works, we plan to further assess the capabilities of the proposed hybrid constraint-handling technique by including other constraints in the asset allocation model. On the one hand, we will consider the so-called risk-budgeting constraints to explicitly control the portfolio risk exposition. On the other hand, based on the European Green Deal and the ESG Disclosure requirements for funds and investments, we will add sustainable-policy constraints to guarantee a minimum level of ESG rating to the investment. In addition, we will consider alternative performance measures in the optimization framework to properly handle tail-risk, and we will extend the experimental part by using different data sets.

A Appendix -Proofs of the main results

Proof of Proposition 3.4: We rewrite problem (3.13) in the following way min

x K : xi∈Ci    1 2 j∈K (x j -y j ) 2    + 1 2 j∈I\K y 2 j ,
that is equivalent to min

x K : xi∈Ci    1 2 j∈K (x j -y j ) 2    - 1 2 j∈K y 2 j + 1 2 j∈I y 2 j . (A.1)
We note that the last term in (A.1) does not depend on x j , j ∈ K, so we can focus our attention on the first two terms, i.e.

min π K (x): xi∈Ci 1 2 π K (x -y) 2 - 1 2 π K (y) 2 .
By contradiction, we suppose that there is a K different from K * , where we recall that K * is the set of indices corresponding to the k largest components of y. At this point, we define 2 . Now, since C i ⊂ R + , we have that ∇f (π K (y)) is non-positive in all components. Moreover, -π K * (y) + π K (y) ≤ 0 due to the fact that π K * (y) is the projection of y onto the set of its largest components. As a result, we obtain g (t) ≥ 0, which implies f (π K (y)) ≥ f (π K * (y)). This means that K * must be the optimal choice.

f (π K (y)) = - 1 2 π K (y) 2 + min π K (x): xi∈Ci
Proof of Proposition 3.5: The proof of the first part of the proposition follows by defining x * K * such that x * K * = δ * ⊗ x * , where ⊗ stands for the Hadamard product. On the contrary, if x * solves (3.16), then x * K * = x * . By taking

δ * i = 1 if i ∈ K * 0 otherwise
we deduce that (δ * , x * ) solves (2.11).

Proof of Proposition 3.7:

The result follows from Theorem 6.27 in [START_REF] Beck | First-order methods in optimization[END_REF], where the hyperplane is represented by the budget constrain (2.5) and the box is [π K * * (l), π K * * (u)]. We recall also that π K * * (z) = (z i1 , . . . , z i k ) with z ij > 0 and i j ∈ K * * .

B Appendix -An approach based on the exact 1 -penalty function

In this appendix we introduce a procedure based on the exact 1 -penalty function for solving cardinalityconstrained portfolio optimization problems. We adapt the approach discussed in [START_REF] Corazza | A novel hybrid PSO-based metaheuristic for costly portfolio selection problems[END_REF] to the algorithms used in the comparison analysis of Subsection 4.1 and thus we define the constraint violations as follows

CV 1 = n i=1 x i -1 CV 2 = max n i=1 δ i -k, 0 CV 3 = n i=1 max {δ i l i -x i , 0} CV 4 = n i=1 max {x i -δ i u i , 0} CV 5 = n i=1 |δ i (1 -δ i )| CV 6 = max n i=1 |x i -x 0,i | -T R, 0 .
In this manner, we introduce the exact 1 -penalty function

F 1 (x, δ; ε) = f (x) + 1 ε 0 [ε 1 CV 1 + ε 2 CV 2 + ε 3 CV 3 + ε 4 CV 4 + ε 5 CV 5 + ε 6 CV 6 ] (B.1)
where ε = (ε 0 , ε 1 , . . . , ε 6 ), with ε > 0 for all i.

The initial parameters vector ε 0 is set to ε 0 = ε 0 0 , ε 0 1 , . . . , ε 0 6 = (10 -4 , 1, . . . , 1) ∈ R 7 , where ε 0 0 is chosen in order to privilege feasible solutions, and the other parameters are equally penalized for all constraint violations.

The vector ε is updated by checking the decrease of the function f (x) and the violation of the constraints. More precisely, on the one hand, every 5 iterations the entry ε 0 (g) is updated according to the rule ε 0 (g + 1) =      min{3 • ε 0 (g), 1} if f (x(g)) ≥ f (x(g -1)) max{0.6 • ε 0 (g), 10 -15 } if f (x(g)) < 0.9 • f (x(g -1)) ε 0 (g) otherwise.

(B.2)

On the other hand, every 10 iterations the entries ε i (g), i = 1, . . . 6, are updated following the scheme

ε i (g + 1) =      min{2 • ε i (g), 10 4 } if CV i (g) > 0.95 • CV i (g -1) max{0.5 • ε g i , 10 -4 } if CV i (g) < 0.9 • CV i (g -1) ε i (g) otherwise, (B.3)
with CV i the respective constraint violation linked to the ε i parameter.

The above quoted strategy privileges optimality of solutions possibly at the expenses of their feasibility, due to the fact that ε 0 (g + 1) in (B.2) is increasing in F 1 (x, δ; ε(g + 1)) when the function value f (x(g)) increases. Moreover, to favour feasibility of solutions possibly at the expenses of their optimality, the penalty parameter ε i (g + 1) in (B.3) is increased when the relative constraint violation in the g-th generation increases with respect to the previous one.

The procedure is also equipped by a splitting and refining technique for the positions of the particles. In particular, at each iteration, a particle p is split in its components x p (g) and δ p (g) that are updated separately. For the vector δ p (g) we employ the following updating rule δ p i (g + 1) = 1 if x p i (g) ∈ [l i , u i ] 0 otherwise for i = 1, . . . , n. Then, δ p (g + 1) is kept fixed and F 1 (x, δ p (g + 1); ε(g + 1)) is minimized with respect to x, obtaining x p (g + 1). Finally, x p (g + 1) is refined getting x p i (g + 1) =

x p i (g + 1)δ p i (g + 1) n i=1 x p i (g + 1)δ p i (g + 1)

, (B.4) for j = 1, . . . , N P .

Lemma 3 . 6

 36 Let l = (l 1 , . . . , l n ) and u = (u 1 , . . . , u n ) be such that l i ≤ u i for i = 1, . . . , n. Let y ∈ R n and define [l, u] = {x ∈ R n : l i ≤ x i ≤ u i }. Then, the orthogonal projection of y onto [l, u] is given component-wise by P [l, u], i (y) = min{max{y i , l i }, u i } (3.17)with i = 1, . . . , n.

Algorithm 2 :

 2 ALLSO-MUT-HInput : x 0 , l, u, k, T R, N L min , N L max , N P , δ, px, ξ Output: x gbest 1 Set g = 0 and N L g = 20 2 Initialize the swarm P g = {x p (g) : p = 1, . . . , N P } and the velocities v p (g) 3 for i = 1 to N P do 4 Project x p (g) onto B using(3.18)

  , 6, 8, 10, 20, 50}, N P = 500, φ = 0.4 [50] ALLSO N L min = 2, N L max = 50, N P = 500, δ = 0.01, px = 0.01, ξ = 10 -6 [44] RLLPSO S = {4, 6, 8, 10, 20, 50}, N P = 500, φ = 0.4, α = 0.4, γ = 0.8, ε = 0.9, ξ = 10 -6 [47] PSO ω min = 0.4, ω max = 0.9, c 1,min = c 2,min = 0.5, c 1,max = c 2,max = 2.5, [37], [15] FA α = 0.5, β min = 0.2, γ = 1 [49]

Figure 1 :

 1 Figure 1: Convergence and diversity analyses on the three data sets. Graphs in the first row show the behaviour of algorithms in terms of mean best value of the objective function, while in the row below are displayed the logarithmic scale plots of the diversity scores.

Figure 2 :

 2 Figure 2: Plots in the first row show the behaviour of the algorithms in terms of mean best value of the penalty function, while in the second row are presented the logarithmic scale graphs of the diversity.

Figure 3 :

 3 Figure 3: Ex-post evolution of net wealth of the benchmark and of the proposed cardinalityconstrained portfolio allocation model with different cardinalities.

1 2 π 1 0g

 121 K (xy)2 andg(t) = f ((1 -t)π K * (y) + tπ K (y)) with t ∈ [0, 1] .Then we have,f (π K (y)) -f (π K * (y)) = g(1) -g(0) = (t) dt and g (t) = ∇f ((1 -t)π K * (y) + tπ K (y)) • (-π K * (y) + π K (y)) ,where ∇f(π K (y)) = -π K (y) + π K (y) -π K (x * ) = -π K (x * ) and π K (x * ) = argmin π K (x): xi∈Ci 1 2 π K (xy)

  .[START_REF] Costa | Theoretical and practical convergence of a self-adaptive penalty algorithm for constrained global optimization[END_REF] The following proposition, whose proof is given in Appendix A, establishes the equivalence between problem (3.16) and the mixed-integer optimization problem(2.11). is a global solution to problem(3.16). Conversely, if x * is a global solution to(3.16), then (δ * , x * ) is a global solution to(2.11), with

Proposition 3.5 We assume that (δ * , x * ) is a global solution to problem (2.11), then x * K *

5

  Calculate the turnover violation using (3.15) Calculate the penalty F for particles in P g 8 Sort P g by F value and divide it in N L g levels 9 Set x gbest (g) = x 1 (g) 10 while g < g max do

	6 end
	7

11 g = g + 1 12

Set P L1 = {x p (g) : p ∈ L 1 } 13 for p = 1 to LP g do

[START_REF] Corazza | Particle swarm optimization with non-smooth penalty reformulation, for a complex portfolio selection problem[END_REF]

Use Algorithm 1 to generate the mutated particle x1, p (g) from x 1, p (g) 15 Project x1, p (g) onto B using (3.18) 16 Calculate the turnover violation using (3.15) 17 end 18 Set P L1 mut = {x p (g) : p ∈ L 1 } 19 Calculate the penalty F for P L1 ∪ P L1 mut and update P L1 based on P L1 mut using F 20 Sort P L1 by the F value 21 Calculate the swarm aggregation indicator using (3.3) and update φ g using (3.7) 22 for p = LP g+1 to N P do 23 Update v p (g) using (3.1) and clamp it using (3.8) 24 Update x p (g) using (3.2) 25 Project x p (g) onto B 26 Calculate the turnover violation using (3.15) 27 end 28 Set Pg be the set of updated particles 29 Calculate φ for P g ∪ Pg and update P g based on Pg using F 30 Sort P g by F value 31 Calculate F for the set x gbest (g), x 1 (g) and update x gbest (g) 32 Calculate the relative improvement t(g) of x gbest (g) using (3.4) 33 if t(g) < t(g -1) or t(g) = 0 then 34 Update N L g using (3.5) and (3.6) 35 end

Table 3 :

 3 Statistics regarding the best values of the objective function over the 25 runs.

	Data set Statistics DLLSO-H DLLSO-MUT-H ALLSO-H ALLSO-MUT-H RLLPSO-H RLLPSO-MUT-H
	S&P 500	mean std	-0.1525 0.0004	-0.1668 0.0015	-0.1524 0.0003	-0.1673 0.0013	-0.1523 0.0004	-0.1624 0.0024

Table 4 :

 4 Relative change of the mutated algorithms versus non-mutated counterparts. The p-values for the paired t-tests are displayed in brackets. Note that in all cases the p-values are under the significance level α = 0.05, indicating the rejection of the null hypothesis of equality of the means, against the alternative left-sided hypothesis.

	Data set	vs.	vs.	vs.
		DLLSO-H (%)	ALLSO-H (%)	RLLPSO-H (%)
	S&P 500	9.3705	9.7948	6.6287
		(6.8333 • 10 -25 )	(4.8584 • 10 -26 )	(7.5173 • 10 -17 )
	Russell 2000	9.3934	10.2666	6.0365
		(2.3999 • 10 -17 )	(2.3049 • 10 -19 )	(7.7679 • 10 -13 )
	Russell 3000	9.2573	10.0354	9.3067
		(2.0685 • 10 -19 )	(1.913 • 10 -23 )	(1.0664 • 10 -17

DLLSO-MUT-H ALLSO-MUT-H RLLPSO-MUT-H

Table 5 :

 5 Comparison with state-of-the-art swarm optimization algorithms implementing the exact 1 -penalty framework proposed in[START_REF] Corazza | A novel hybrid PSO-based metaheuristic for costly portfolio selection problems[END_REF].

	Data set Statistics	PSO-1	FA-1	ALLSO-MUT-1 ALLSO-MUT-H
	S&P 500	feasible sol. (%) mean CV mean F 1	0 3.1572 • 10 -12 3.1572 • 10 -12 0 -0.1325 -0.1325	100 0 -0.1325	100 0 -0.1673
	Russell 2000	feasible sol. (%) mean CV mean F 1	68 1.6035 • 10 -5 -0.1718	100 0 -0.1639	100 0 -0.1639	100 0 -0.2122
	Russell 3000	feasible sol. (%) mean CV mean F 1	80 7.4015 • 10 -12 -0.1866	100 0 -0.1291	0 0.0320 -0.1493	100 0 -0.2300

Table 6 :

 6 Structure of transaction costs.

Table 7 :

 7 Performance of the proposed cardinality-constrained portfolio allocation model for different cardinalities in comparison to the benchmark.

	k %	30%	15%	5%	2%	Benchmark
	num. assets	335	167	55	22	1119
	Sharpe ratio 0.4917	0.4794	0.4467	0.2878	0.1731
	Omega ratio 3.3269	3.2428	3.0360	2.0980	1.5388
	CAGR	1.1148	1.6899	2.1118	2.2502	1.3460
	std	0.0052	0.0052	0.0054	0.0077	0.0066
	mean DD	-0.0066 -0.0057 -0.0050 -0.0072	-0.0070
	std DD	0.0070	0.0064	0.0060	0.0106	0.0089
	mean λ ($)	16,046	10,713 6,571.7 3,121.9	-
	Λ %	0.1536	0.1011	0.0610	0.0290	-

We denote by x the floor of x and by x%y the rest of the division of x by y.

The gross return of asset i at month t is defined as R g i,t = S i,t S i,t-1 , where Si,t is the price of the i-th asset at the end of month t.

Acknowledgements

The third author is member of the INdAM (Italian Institute for Advanced Mathematics) group.