domain: phys stat

text:
Amelia Carolina Sparavigna

Relations Between Tsallis and Kaniadakis Entropic Measures and Rigorous Discussion of Conditional Kaniadakis Entropy

Keywords:

Relations Between Tsallis and Kaniadakis Entropic Measures and Rigorous Discussion of Conditional

Kaniadakis Entropy

Introduction

In information theory, measures of information can be obtained from the probability distribution of some events contained in a sample set of possible events. These measures are the entropies.

In 1948 [START_REF] Shannon | A Mathematical Theory of Communication[END_REF], Claude Shannon defined the entropy H (Greek letter Eta), of a discrete random variable X, as the expected value of the information content: H(X)

∑ = i i i ) )I(x p(x ∑ - = i i b i ) p(x log ) p(x [2].
In this expression, I is the information content of X, the probability of i-event is i p and b is the base of the used logarithm. Common values of the base are 2, Euler's number e, and 10. Besides Shannon entropy, several other entropies are used in information theory; here we will consider the generalized entropies of Tsallis and Kaniadakis (also known as κ-entropy) entropies [START_REF] Borda | Fundamentals in Information Theory and Coding[END_REF][START_REF] Tsallis | Possible Generalization of Boltzmann-Gibbs Statistics[END_REF][START_REF] Kaniadakis | Statistical Mechanics in the Context of Special Relativity[END_REF]. We will show relations between them. We will also propose a rigorous discussion of the conditional Kaniadakis entropy, deduced from these relations.

The entropies

In the following formulas we can see how the abovementioned entropies (Shannon, Tsallis and Kaniadakis) are defined, with a corresponding choice of measurement units equal to 1:

i i i p ln p S : Shannon (1) ∑ - =         ∑ - - = = i q i q p 1 1 q 1 T T : Tsallis (2) ∑ - - = - + i κ 1 i κ 1 i κ 2κ p p K : entropy) - (κ Kaniadakis (3)
In ( 2),(3) we have entropic indices q and κ , and:

S K 0 κ lim ; S T 1 q lim = → = → .
Let us consider the joint entropy H(A,B) of two independent subsystems A,B. We have:

S(B) S(A) B) S(A, (4)

+ = q)T(A)T(B) (1 T(B) T(A) B) T(A, (5) - + + = ∑ + = ℑ ℑ + ℑ = - + i κ 1 i κ 1 i 2 p p with (A) K(B) (B) K(A) B) K(A, (6)
Note that for the generalized additivity of κ -entropy, we need another function containing probabilities (see [START_REF] Sparavigna | On the Generalized Additivity of Kaniadakis Entropy[END_REF] and references therein).

Basic relations between K and T

Let us consider κ -entropy K. We have that:

DOI: 10.18483/ijSci.866 , 2 T T K (7) κ 1 κ 1 κ - + + =
In (7) we used the Tsallis entropies:

κ 1 p κ 1 κ) 1 T(q ; κ 1 p κ 1 κ) 1 T(q i κ 1 i i κ 1 i - ∑ = - = + ∑ - = + = - +
In fact, we have:

            ∑ - ∑ - = - ∑ + + ∑ - = - + - + i κ 1 i i κ 1 i i κ 1 i i κ 1 i p p 2κ 1 κ 1 p κ 1 κ 1 p κ 1 2 1 K
Eq.( 7) is a simpler form of an expression given in [START_REF] Kaniadakis | Non-linear kinetics underlying generalized statistics[END_REF][START_REF] Santos | Generalized quantum entropies[END_REF]. However, besides this relation, because of the generalized additivity possessed by the Kaniadakis entropy, we need also another relation, concerning function ℑ . It is:

        + + - = ℑ - + κ 2 T T 2 κ (8) κ 1 κ 1 κ
In fact:

∑ + = + - ∑ + - ∑ = ℑ - + - +       i κ 1 i κ 1 i i κ 1 i i κ 1 i 2 p p κ 2 κ 1 p κ 1 κ 1 p κ 1 2
κ In ( 7) and ( 8), we have the Kaniadakis functions expressed by the Tsallis entropy. However, we can also write T expressed by means of Kaniadakis entropy.

        + + + - + = ℑ + - + - + κ 2 T T κ 2 2K κ 1 κ 1 κ 1 κ 1 T T
And then:

κ 1 T K κ 1 κ κ κ 1 (9) + - = ℑ + Let us have: q 1 κ - = . q) (1 1 K (10) q 1 q 1 q T - - ℑ + - - =
We can have also:

        + + + - - = ℑ - - + - + κ 2 T T κ 2 2K κ 1 κ 1 κ 1 κ 1 T T So that: κ 1 T K κ 1 κ κ κ 1 (11) - + = ℑ - Let us have: 1 κ - = q
. We have again Eq.10.

Generalized additivity

Let us consider two subsystems A and B. We can find a relation between the joint Tsallis and Kaniadakis entropies. Using [START_REF] Abe | Nonadditive Conditional Entropy and its Significance for Local Realism[END_REF], we obtain:

q 1 1 q 1 B) (A, B) (A, K B) (A, T (12)

  When the subsystems are independent, for Tsallis entropy we have Eq.5, and then: DOI: 10.18483/ijSci.866

To continue, we can assume:

This relation was proposed in [START_REF] Scarfone | Thermodynamic equilibrium and its stability for microcanonical systems described by the Sharma-Taneja-Mittal entropy[END_REF], but it is clear that it can obtained from (13):

, we must have (14) to fulfil (13).

Conditional entropy

In a previous paper [START_REF] Sparavigna | Mutual Information and Nonadditive Entropies: A Method for Kaniadakis Entropy[END_REF] we proposed, in the framework of a discussion on mutual entropy, an expression for the conditional Kaniadakis entropy. If the entropic index has a low value, the formula we find in [START_REF] Sparavigna | Mutual Information and Nonadditive Entropies: A Method for Kaniadakis Entropy[END_REF] can be considered an approximation of the expression that we are here deducing. Let us start from the conditional Tsallis entropy [START_REF] Abe | Nonadditive Conditional Entropy and its Significance for Local Realism[END_REF]:

. From now on, we do not write the indices of Kaniadakis functions. However, remember that they exist and have value (1-q). Let us use [START_REF] Abe | Nonadditive Conditional Entropy and its Significance for Local Realism[END_REF], written for joint and conditional entropies:

We have:

From this equation we find:

We can divide (18) in two equations:

Let us note that ( 19) and ( 20) are generalizing ( 6) and ( 14). In both ( 19),(20) we have conditional functions

. Then, the rigorous expression of Kaniadakis conditional entropy is:

In the case we are using the entropic index with a low value, this expression is approximated by:

The discussion in [START_REF] Sparavigna | Mutual Information and Nonadditive Entropies: A Method for Kaniadakis Entropy[END_REF] can be proposed again using (21). For instance, if we consider the mutual entropy (without renormalization): and then the mutual information is zero. Moreover, when we have a small value of the entropic index, function ℑ is equal to 1.


R Pedro Brito

Luís N Vicente

Efficient cardinality/mean-variance portfolios

Keywords: portfolio selection, cardinality, sparse portfolios, multiobjective optimization, efficient frontier, derivative-free optimization

We propose a novel approach to handle cardinality in portfolio selection, by means of a biobjective cardinality/mean-variance problem, allowing the investor to analyze the efficient tradeoff between returnrisk and number of active positions. Recent progress in multiobjective optimization without derivatives allow us to robustly compute (in-sample) the whole cardinality/mean-variance efficient frontier, for a variety of data sets and mean-variance models. Our results show that a significant number of efficient cardinality/mean-variance portfolios can overcome (out-of-sample) the naive strategy, while keeping transaction costs relatively low.

Introduction

One knows since the pioneer work of Markowitz [START_REF] Markowitz | Portfolio selection[END_REF] that a rational investor has typically two goals in mind: to maximize the portfolio return (given, e.g., by the portfolio expected return) and to minimize the portfolio risk (described, e.g., by the portfolio variance). Traditionally, the Markowitz mean-variance optimization model is taken as a quadratic program (QP), intended to minimize the portfolio risk (variance) for a given level of expected return, over a set of feasible portfolios. By varying the level of expected return, the Markowitz model determines the so-called efficient frontier, as the set of nondominated portfolios regarding the two goals (variance and mean of the return). The rational investor can thus make choices, by analyzing the tradeoff between expected return and variability of the investment, over a set of appropriate portfolios.

Several modifications to the classical Markowitz model or alternative methodologies have since then been proposed. One resulting from a simple observation was suggested in an article by DeMiguel, Garlappi, and Uppal [START_REF] Demiguel | Optimal versus naive diversification: How inefficient is the 1/N portfolio strategy?[END_REF]. These authors analyzed a number of methodologies inspired on the classic model of Markowitz and showed that none were able to significantly and consistently overcome the naive strategy, that is to say, the one in which the available investor's wealth is divided equally among the available securities. One possible explanation is related to the ill conditioning of the objective function of the Markowitz model (given by the variance of the return).

One of the important issues to consider in portfolio selection is how to handle transaction costs. There are well known modifications that can be made in the Markowitz model to incorporate transaction costs, such as to bound the turnover, which basically amount to further linear constraints in the QP. A recent technique to keep transaction costs low consists of selecting sparse portfolios, i.e., portfolios with few active positions, by imposing a cardinality constraint. Such a constraint, however, changes the classical QP into a MIQP (mixed-integer quadratic programming), which can no longer be solved in polynomial time.

In this paper, we suggest an alternative approach to the cardinality constrained Markowitz mean-variance optimization model, reformulating it directly as a biobjective problem, allowing the investor to analyze the tradeoff between cardinality and mean-variance, in a general scenario where short-selling is permitted. Such an approach allows us to find the set of nondominated points of biobjective problems in which an objective is smooth and combines mean and variance and the other is nonsmooth (the cardinality or 0 norm of the vector of portfolio positions). The mean-variance objective function can take a number of forms. A parameter free possibility is given by profit per unity of risk (a nonlinear function obtained by dividing the expected return by its variance).

Given the lack of derivatives of the cardinality function, we decided then to apply a directional derivative-free algorithm for the solution of the biobjective optimization problem. Such methods do not require derivatives, although their convergence results typically assume some weak form of smoothness such as Lipschitz continuity. Direct multisearch is a derivative-free multiobjective methodology for which one can show some type of convergence in the discontinuous case. More importantly, it exhibited excellent numerical performance on a comparison to a number of other multiobjective optimization solvers. We applied direct multisearch to determine (in-sample) the set of efficient or nondominated cardinality/mean-variance portfolios.

To illustrate our approach, we gathered several data sets from the FTSE 100 index (for returns of single securities) and from the Fama/French benchmark collection (for returns of portfolios), computed the efficient cardinality/meanvariance portfolios using (in-sample) optimization, and measured their out-ofsample performance using a rolling-sample approach. We found that a large number of sparse portfolios for the FTSE 100 data sets, among the efficient cardinality/mean-variance ones, consistently overcome the naive strategy in terms of out-of-sample performance measured by the Sharpe ratio. This effect is also clearly visible for the FF data sets, where the performance of a large portion of the cardinality/mean-variance efficient frontier outperforms, in most of the instances, the naive strategy. The transactions costs are shown to be relatively low for all efficient cardinality/mean-variance portfolios, with a moderate increase with cardinality.

The organization of our paper is as follows. In the next section, we formulate the classical Markowitz model for portfolio selection, describe the naive strategy, and formulate the problem with cardinality constraint. In Section 3, we reformulate the cardinality constrained Markowitz mean-variance optimization model as a biobjective problem for application of multiobjective optimization. In Section 4, we present the empirical results. Finally, in Section 5 we summarize our findings and discuss future research.

Portfolio selection models

The classical Markowitz mean-variance model

Portfolios consist of securities (shares or bonds, for example, or classes or indices of the same). Suppose the investor has a certain wealth to invest in a set of N securities. The return of each security i is described by a random variable R i , whose average can be computed (from estimation based on historical data). Let µ i = E(R i ), i = 1, . . . , N , denote the expected returns of the securities. Let also w i , i = 1, . . . , N , represent the proportions of the total investment to allocate in the individual securities. The portfolio return is assumed linear in w 1 , . . . , w N , and thus the portfolio expected return can be written as

E(R) = E(w i R 1 + • • • + w N R N ) = w 1 µ 1 + • • • + w N µ N = µ w with µ = (µ 1 , . . . , µ N )
and w = (w 1 , . . . , w N ) .

The portfolio variance, in turn, is calculated by

V (R) = E N i=1 w i R i -E N i=1 w i R i 2 .
So,

V (R) = N i=1 N j=1 E[(R i -µ i )(R j -µ j )]w i w j .
Representing each entry i, j of the covariance matrix Q by

σ ij = E[(R i -µ i )(R j -µ j )],
one has

V (R) = w Qw,
where Q is symmetric and positive semi-definite (and typically assumed positive definite). As said before, a portfolio is defined by an N × 1 vector w of weights representing the proportion of the total funds invested in the N securities. This vector of weights is thus required to satisfy the constraint

N i=1 w i = e w = 1,
where e is the N × 1 vector of entries equal to 1. Lower bounds on the variables, of the form w i ≥ 0, i = 1, . . . , n, can be also considered if short selling is undesirable. In general, we will say that L i ≤ w i ≤ U i , i = 1, . . . , N , for given lower L i and upper U i bounds on the variables. Markowitz's model [START_REF] Markowitz | Portfolio selection[END_REF][START_REF] Markowitz | Portfolio selection: Efficient diversification of investments[END_REF] is based on the formulation of a mean-variance optimization problem. By solving this problem, we identify a portfolio of minimum variance among all which provide an expected return not below a certain target value r. The aim is thus to minimize the risk from a given level of return. The formulation of this problem can be described as:

min w∈R N w Qw subject to µ w ≥ r, e w = 1, L i ≤ w i ≤ U i , i = 1, . . . , N. (1)
Problem ( 1) is a convex quadratic programming problem (QP), for which the first order necessary conditions are also sufficient for (global) optimality. See [START_REF] Cornnuejols | Optimizations Methods in Finance[END_REF][START_REF] Steinbach | Markowitz revisited: Mean-variance models in financial portfolio analysis[END_REF] for a survey of portfolio optimization. The classical Markowitz mean-variance model can be seen as way of solving the biobjective problem which consists of simultaneously minimizing the portfolio risk (variance) and maximizing the portfolio profit (expected return)

L i ≤ w i ≤ U i , i = 1, . . . , N. (2)
In fact, it is easy to prove that a solution of (1) is nondominated, efficient or Pareto optimal for (2). Efficient portfolios are thus the ones which have the minimum variance among all that provide at least a certain expected return, or, alternatively, those that have the maximal expected return among all up to a certain variance. The efficient frontier (or Pareto front) is typically represented as a 2-dimensional curve, where the axes correspond to the expected return and the standard deviation of the return of an efficient portfolio.

The naive strategy 1/N

The naive strategy is the one in which the available investor's wealth is divided equally among the securities available

w i = 1 N , i = 1, . . . , N.
This strategy has diversification as its main goal, it does not involve optimization, and it completely ignores the data. Although a number of theoretical models have been developed in the last years, many investors pursuing diversification revert to the use of the naive strategy to allocate their wealth (see [START_REF] Benartzi | Naive diversification strategies in defined contribution saving plans[END_REF]). DeMiguel, Garlappi, and Uppal [START_REF] Demiguel | Optimal versus naive diversification: How inefficient is the 1/N portfolio strategy?[END_REF] evaluated fourteen models across seven empirical data sets and showed that none is consistently better than the naive strategy. A possible explanation for this phenomenon lies on the fact that the naive strategy does not involve estimation and promotes 'optimal' diversification. The naive strategy is therefore an excellent benchmarking strategy.

The cardinality constrained Markowitz mean-variance model

Since the appearance of the classical Markowitz mean-variance model, a number of methodologies have been proposed to render it more realistic. The classical Markowitz model assumes a perfect market without transaction costs or taxes, but such costs are an important issue to consider as far as the portfolio selection is concerned, especially for small investors. Recently, it has been studied the addition of a constraint that sets an upper bound on the number of active positions taken in the portfolio, in an attempt to improve performance and reduce transactions costs. Such a cardinality constraint is defined by limiting card(x) = |{i ∈ {1, ..., N } : x i = 0}| and leads to cardinality constrained portfolio selection problems. In particular, the cardinality constrained Markowitz mean-variance optimization problem has the form:

L i ≤ w i ≤ U i , i = 1, . . . , N, (3)
where K ∈ {1, . . . , N }. Although card(x) is not a norm, it is frequently called the 0 norm in the literature, x 0 = card(x). By introducing binary variables, one can rewrite the problem as a mixed-integer quadratic programming (MIQP) problem: min

w,y∈R N w Qw subject to µ w ≥ r, e y ≤ K, e w = 1, L i y i ≤ w i ≤ U i y i , i = 1, . . . , N, y i ∈ {0, 1}, i = 1, . . . , N. (4)
However such MIQPs are known to be hard combinatorial problems. The number of sparsity patterns in w (i.e., number of different possibilities of having K

nonzeros entries) is N K = N !/[(N -K)!K!].
Although there are exact algorithms for the solution of MIQPs (see [START_REF] Bertsimas | Algorithm for cardinality-constrained quadratic optimization[END_REF][START_REF] Bienstock | Computational study of a family of mixed-integer quadratic programming problems[END_REF][START_REF] Boyd | Convex Optimization[END_REF][START_REF] Vielma | A lifted linear programming branchand-bound algorithm for mixed-integer conic quadratic programs[END_REF]), many researchers and portfolio managers prefer to use heuristics approaches (see [START_REF] Bach | Convex relaxations for subset selection[END_REF][START_REF] Cesarone | Efficient algorithms for mean-variance portfolio optimization with hard real-world constraints[END_REF][START_REF] Chang | Heuristics for cardinality constrained portfolio optimisation[END_REF][START_REF] Fieldsend | Cardinality constrained portfolio optimisation[END_REF][START_REF] Lin | A multiobjective genetic algorithm for portfolio selection[END_REF][START_REF] Woodside-Oriakhi | Heuristic algorithms for the cardinality constrained efficient frontier[END_REF]). Some of these heuristics vary among evolutionary algorithms, tabu search, and simulated annealing (see [START_REF] Fieldsend | Cardinality constrained portfolio optimisation[END_REF][START_REF] Woodside-Oriakhi | Heuristic algorithms for the cardinality constrained efficient frontier[END_REF]).

Promotion of sparsity is also used in the field of signal and imaging processing, where a new technique called compressed sensing has been intensively studied in the recent years. Essentially one aims at recovering a desired signal or image with the least possible amount of basis components. The major developments in compressed sensing have been achieved by replacing the 0 norm by the 1 one, the latter being a convex relation of the former and known to also promote sparsity. The use of the 1 norm leads to recovering optimization problems solvable in polynomial time (in most of the cases equivalent to linear programs), and a number of sparse optimization techniques have been developed for the numerical solution of such problems. These ideas have already been used in portfolio selection primarily to promote regularization of ill conditioning (of the estimation of data or of the variance of the return itself). DeMiguel et al. [START_REF] Demiguel | A generalized approach to portfolio optimization: Improving performance by constrained portfolio norms[END_REF] constrained the Markowitz classical model by imposing a bound on the 1 norm of the vector of portfolio positions, among other possibilities. Brodie et al. [START_REF] Brodie | Sparse and stable Markowitz portfolios[END_REF] focus on a modification to the Markowitz mean-variance classical model by the incorporation of a term involving a multiple of the 1 norm of the vector of portfolio positions. Inspired by sparse reconstruction (see, for instance, [START_REF] Boyd | Convex Optimization[END_REF]), they also proposed an heuristic for the solution of the problem.

The cardinality/mean-variance biobjective model

Although the cardinality constrained Markowitz mean-variance model described in [START_REF] Bach | Convex relaxations for subset selection[END_REF] provides an alternative to the classical Markowitz model in the sense of realistically limiting the number of active positions in a portfolio, it is dependent on the parameter K, the maximum number of such positions. Thus, one has to vary K to obtain various levels of cardinality or sparsity, and for each value of K solve an MIQP of the form (4).

The alternative suggested in this paper is to consider the cardinality function as an objective function itself. At a first glance, one could see the problem as a triobjective optimization problem by minimizing the variance of the return, maximizing the expected return, and minimizing the cardinality over the set of feasible portfolios. Such a framework was taken into account in the studies [START_REF] Anagnostopoulos | A portfolio optimization model with three objectives and discrete variables[END_REF][START_REF] Anagnostopoulos | The mean-variance cardinality constrained portfolio optimization problem: An experimental evaluation of five multiobjective evolutionary algorithms[END_REF][START_REF] Cesarone | A new method for mean-variance portfolio optimization with cardinality constraints[END_REF][START_REF] Lorenzo | A concave optimization-based approach for sparse portfolio selection[END_REF]. However, these authors did not investigate the effects of cardinality constraints on portfolio models in terms of out-of-sample performance, a subject still poorly analyzed in the literature. On the other hand, investors may find it useful to directly analyze the tradeoff between cardinality and meanvariance. A parameter-free possibility is to consider a Sharpe ratio type objective function, by maximizing expected return per variance and minimizing the cardinality, over the set of feasible portfolios. In this case, the cardinality/mean-variance biobjective optimization problem is posed as

min w∈R N -µ w w Qw min w∈R N card(w)
subject to e w = 1,

L i ≤ w i ≤ U i , i = 1, . . . , N. (5)
By solving [START_REF] Bertsimas | Algorithm for cardinality-constrained quadratic optimization[END_REF], we identify a cardinality/mean-variance efficient frontier. A portfolio in this frontier is such that there exists no other feasible one which simultaneously presents a lower cardinality and a lower mean-variance measure. Given such an efficient frontier and a mean-variance target, an investor may directly find the answers to the questions of what is the optimal (lowest) cardinality level that can be chosen and what are the portfolios leading to such a cardinality level. Problem (5) has two objective functions and linear constraints. The first objective f 1 (w) = -µ w/w Qw is nonlinear but smooth. However, the second objective function f 2 (w) = card(w) = |{i ∈ {1, ..., N } : w i = 0}| is piecewise linear discontinuous, consequently nonlinear and nonsmooth. We have thus decided to solve the biobjective optimization problem (5) using a derivative-free solver, based on direct multisearch.

Empirical performance of efficient cardinality/meanvariance portfolios

Now we report a number of experiments made to numerically determine and assess the efficient cardinality/mean-variance frontier. We applied direct multisearch to determine the Pareto front or efficient frontier of the biobjective optimization problem (5) (according to Appendix A). We tested three data sets collected from the FTSE 100 index and three others from the Fama/French benchmark collection (see Subsection 4.1). The efficient frontiers obtained by the initial in-sample optimization are given in Subsection 4.2.

The out-of-sample performance of the cardinality/mean-variance efficient portfolios, measured by a rolling-sample approach, is described in Subsection 4.3. In Subsection 4.4 we measure the out-of-sample performance by the Sharpe ratio, in Subsubsection 4.5 we report the proportional transaction costs, and in Subsection 4.6 we measure the out-of-sample performance by the Sharpe ratio of returns net of transaction costs, all of this for each cardinality/mean-variance efficient portfolio. To better assess the robustness of our results, we also considered, using the FTSE 100 data, a sample including the financial crisis years 2008-2010, and the corresponding results are reported in Subsection 4.7. The section is ended with a discussion of the overall obtained results.

Data sets

For the first three data sets we collected daily data for securities from the FTSE 100 index, from 01/2003 to 12/2007 (five years). Such data is public and available from the site http://www.bolsapt.com. The three data sets are referred to as DTS1, DTS2, and DTS3, and are formed by 12, 24, and 48 securities, respectively. The composition of these data sets is given in Table 1. We used the daily continuous returns for the in-sample optimization (estimation of Q and µ) and the daily discrete returns for the out-of-sample analysis. We also included in our experiments three data sets from the Fama/French benchmark collection (FF10, FF17, and FF48, with cardinalities 10, 17, and 48), using the monthly returns from 07/1971 to 06/2011 (forty years) given there for a number of industry security sectors. More information on these security sectors (or portfolios of securities) can be found in http://mba.tuck.dartmouth.edu/pages/faculty/ ken.french/data_library.html.

SECURITIES 3 I GROUP (1,2,3) JOHNSON MATTHEY P (3) AMEC (1,2,3) LEGAL & GENERAL (3) ANGLO AMERICAN (1,2,3) LLOIDS BANKING GR (3) ANTOFAGASTA (1,2,3) LONMIN (3) ASSOCIAT BRIT FOO (1,2,3) MARKS & SPENCER (3) ASTRAZENECA (1,2,3) MORRINSON SUPERMKT (3) AVIVA (1,2,3) NEXT (3) B SKY B GROUP (1,2,3) OLD MUTUAL (3) BAE SYSTEMS (1,2,3) PEARSON (3) BARCLAYS (1,2,3) PRUDENTIAL (3) BG GROUP (1,2,3) REED ELSEVIER PLC (3) BHP BILLITON (1,2,3) RENTOKIL INITIAL (3) BP (2,3) REXAM (3) BRIT AMER TOBACCO (2,3) RIO TINTO (3) BRIT LAND CO REIT (2,3) ROYAL BK SCOTL GR (3) BRITISH AIRWAYS (2,3) RSA INSUR GRP (3) CAB & WIRE WRLD (2,3) SABMILLER (3) CAPITA GRP (2,3) SAGE GRP (3) COBHAM (2,3) SAINSBURY (3) DIAGEO (2,3) SCHRODERS (3) HAMMERSON REIT (2,3) SEVERN TRENT (3) IMPERIAL TOBACCO (2,3) SHIRE (3) INTERNATIONAL POW (2,3) UNITED UTILITIES (3) INVENSYS (2,3) VODAFONE GRP (3)
Table 1. Composition of the three data sets from the FTSE 100 index. In brackets we indicate the data set to which each security belongs to.

In-sample optimization

We then applied the solver dms (version 0.2) to compute the efficient frontier (or Pareto front) of the cardinality/mean-variance biobjective optimization problem [START_REF] Bertsimas | Algorithm for cardinality-constrained quadratic optimization[END_REF]. A few modifications to (5) were made before applying the solver as well as a few changes to the solver default parameters (the details are described in Appendix A). We present results for the initial in-sample optimization.

This instance was solved using the quadprog function from the MATLAB [START_REF] Mathworks | [END_REF] Optimization Toolbox. Regarding problem [START_REF] Boyd | Convex Optimization[END_REF], it is known that not allowing short-sale has a regularizing effect on minimum-variance Markowitz portfolio selection (see [START_REF] Jagannathan | Risk reduction in large portfolios: Why imposing the wrong constraints hfelps[END_REF]) and leads to portfolios of low cardinality. Since we know that minimum variance portfolios outperform mean-variance portfolios (the estimate error of the expected returns is eliminated, see [START_REF] Jagannathan | Risk reduction in large portfolios: Why imposing the wrong constraints hfelps[END_REF]), we considered the following cardinality constrained minimum variance model (instead of the one introduced in Section 2.3)

min w∈R N w Qw subject to card(w) ≤ K, e w = 1, L i ≤ w i ≤ U i , i = 1, . . . , N.
L i y i ≤ w i ≤ U i y i , i = 1, . . . , N, y i ∈ {0, 1}, i = 1, . . . , N. (8)
We also mark in the plots the portfolios that result from solving problem [START_REF] Brodie | Sparse and stable Markowitz portfolios[END_REF] for each value of K ∈ [1, N ]. For this purpose we used the solver cplexmiqp from ILOG IBM CPLEX for MATLAB [START_REF] Ibm | IBM ILOG CPLEX[END_REF].

Out-of-sample performance

The analysis of out-of-sample performance relies on a rolling-sample approach. For the FTSE 100 data sets we considered 12 periods (months) of evaluation. We begin by computing the efficient frontier (or Pareto front) of the cardinality/mean-variance biobjective optimization problem (5) for the in-sample time window from 01/2003 to 12/2006 (see Subsection 4.2). We then held fixed each portfolio and observed its returns over the next period (January 2007). Then we discarded January 2003 and brought January 2007 into the sample. We repeated this process until exhausting the 12 months of 2007. We applied the same rolling-sample approach to the FF data sets, considering an initial in-sample time window from 07/1971 to 06/1996 (see Subsection 4.2) and 15 periods of evaluation (the 15 next years).

Out-of-sample performance measured by the Sharpe ratio

In each period of evaluation, the out-of-sample performance was then measured by the Sharpe ratio

S = m -r f σ ,
where m is the mean return, r f is the return of the risk-free asset3 , and σ is the standard deviation. The results (over all the periods of evaluation) are given in Figures 789for the FTSE 100 portfolios and in Figures 101112for the FF ones. Using IBM SPSS Statistics [23] we calculated the p-values for the statistical significance of the difference between Sharpe ratios of the benchmark naive portfolio and all the others computed portfolios. We did not report them here because they are not statistically significant.

Transaction costs

Since one is rebalancing portfolios for each out-of-sample period, one can compute the transaction costs of such a trade. We set the proportional transaction cost equal to 50 basis points per transaction (as usually assumed in the literature). Thus the cost of a trade over all assets is given by

T C = T -1 t=1 0.5% N i=1 | w i,t+1 -w i,t |, (9)
with T = 12 for the FTSE 100 data sets and T = 15 for the FF data sets. The results are given in Figures 131415for the FTSE 100 portfolios and in Figures 16-18 for the FF ones.

Out-of-sample performance measured by the Sharpe ratio of returns net of transaction costs

In the presence of transaction costs we calculated the Sharpe ratio of returns net of transaction costs

SR = m -T C -r f σ ,
where m is the mean return, T C is the proportional transaction cost in (9), r f is the return of the risk-free asset, and σ is the standard deviation. The outof-sample performance was then measured by the Sharpe ratio of returns net of transaction costs. The results are given in ). The data sets were formed as described in Section 4.1, but excluding British Airways (see Table 1) due to missing data during the period considered, and including Wolseley (following an arbitrary alphabetic order). We performed an out-of-sample analysis as described in Section 4.3. We used daily periods of evaluation. We began by computing the efficient frontier (or Pareto front) of the cardinality/mean-variance biobjective optimization problem (5) for the insample time window from 01/2003 to 12/2010 (using daily data). We then held fixed each portfolio and observed its returns over the next period (first trading day of January 2011). Then we discarded this first trading day of January 2011 and brought this into the sample. We repeated this process until exhausting the firsts 15 trading days of 2011.

The results of the out-of-sample performance measured by the Sharpe ratio4 are given in Figures 252627. The results of the proportional transaction costs, are given in Figures 282930. The results of the out-of-sample performance measured by the Sharpe ratio of returns net of transaction costs, are given in Figures 313233.

Discussion of the results

Contrary to one could think, given the intractability of f 2 (w) = card(w) and the fact that no derivatives are being used for f 1 (w) = -µ w/w Qw, direct multisearch (the solver dms) was capable of quickly determining (in-sample) the efficient frontier for the biobjective optimization problem [START_REF] Bertsimas | Algorithm for cardinality-constrained quadratic optimization[END_REF]. For instance, for the data sets of roughly 50 assets, a regular laptop takes a few dozens of seconds to produce the efficient frontiers. We have a direct way of dealing with sparsity, which offers a complete determination of an efficient frontier for all cardinalities. According to a priori preferences, one could choose (in-sample) the desired cardinality. For the portfolios constructed using the FTSE 100 index data (portfolios of individual securities), a large number of our sparse portfolios, among the efficient cardinality/mean-variance ones, consistently overcame the naive strategy and at least one of the two related classical Markowitz models, in terms of out-of-sample performance measured by the Sharpe ratio. This effect has even happened for the largest data set (DTS3 with 48 securities), where the demand for sparsity is more relevant. For the portfolios constructed using the Fama/French benchmark collection (where securities are portfolios rather than individual securities), the scenario is different since the behavior of the naive strategy is even more difficult to outperform. Still, a large number of sparse efficient cardinality/mean-variance portfolios consistently overcame the naive strategy.

In both cases, FTSE 100 and FF data, the transaction costs of the efficient cardinality/mean-variance portfolios are lower than the mean per variance portfolio (solution of problem ( 6)) and higher than the minimum-variance portfolio (solution of problem ( 7)). Note that the minimum-variance portfolio does not allow short-selling, and so the weights at the outset are much more limited, thus leading to better results. Evaluating the performance out-of-sample by the Sharpe ratio of returns net of transaction costs (take into account the transaction costs), the efficient cardinality/mean-variance portfolios do not overcame the naive strategy for FTSE 100 data, but for FF data a large number of sparse efficient cardinality/mean-variance portfolios still consistently overcame the naive strategy. When we compare the performance results between the efficient cardinality/mean-variance portfolios and the cardinality constrained minimum variance portfolios (solution of ( 8)), without considering the transaction costs, we observed better results for the FTSE 100 and worse for the FF. The MIQP performed better in terms of Sharpe ratio of returns net of transaction costs since the cost of transaction costs are lower, one possible explanation for this is the fact of not taking into account the estimation of the expected returns. Moreover, our cardinality/mean-variance portfolios are truly efficient whereas the cardinality constrained minimum variance do not necessarily exhibit Pareto efficiency. For the FTSE 100 data set, the analysis including the financial crisis years 2008-2010 shows that the results are robust. Finally, we also computed the cardinality/mean-variance efficient frontier for the data set FF100, where portfolios are formed on size and book-to-market (see Figure 34). (This time we needed a budget of the order of 10 7 function evaluations, see Appendix A.) We remark that FF48 and FF100 are the data sets also used in [START_REF] Brodie | Sparse and stable Markowitz portfolios[END_REF]. In this paper, as we said before, the authors focus on a modification to the Markowitz classical model by the incorporation of a term involving a multiple of the 1 norm of the vector of portfolio positions. Despite the different sparse-oriented techniques and different strategies for evaluating out-of-sample performance, in both approaches (theirs and ours), sparse portfolios are found overcoming the naive strategy. In our approach one computes sparse portfolios satisfying an efficient or nondominant property and one does it directly and in single run, whereas in [START_REF] Brodie | Sparse and stable Markowitz portfolios[END_REF], there is a need to vary a tunable parameter and select the portfolios according to some criterion to be met (for example, sparsity). It is unclear what sort of efficient or nondominant property their portfolios satisfy. Moreover, we provide results for all cardinality values (from 1 to 48 in FF48 and from 1 to 100 in FF100), while in [START_REF] Brodie | Sparse and stable Markowitz portfolios[END_REF] the authors report results for cardinality values from 4 and 48 (FF48) and from 3 to 60 (FF100). We therefore claim to have a more direct way of dealing with sparsity, which offers a complete determination of an efficient frontier for all cardinalities.

Conclusions and perspectives for future work

In this paper we have developed a new methodology to deal with the computation of mean-variance Markowitz portfolios with pre-specified cardinalities. Instead of imposing a bound on the maximum cardinality or including a penalization or regularization term into the objective function (in classical Markowitz mean-variance models), we took the more direct approach of explicitly considering the cardinality as a separate goal. This led us to a cardinality/mean-variance biobjective optimization problem (5) whose solution is given in the form of an efficient frontier or Pareto front, thus allowing the investor to tradeoff among these two goals when having transaction costs and portfolio management in mind. In addition, and surprisingly, a significant portion of the efficient cardinality/meanvariance portfolios (with cardinality values considerably lower than the number N of securities) have exhibited superior out-of-sample performance (under reasonably low transaction costs that only increase moderately with cardinality). We solved the biobjective optimization problem (5) using a derivative-free solver running direct multisearch. Direct-search methods based on polling are known in general to be slow but extremely robust due their directional properties. Such a feature is crucial given the difficulty of the problem (one discontinuous objective function, the cardinality, and discontinuous Pareto fronts). We have observed the robustness of direct multisearch, in other words, its capability of successfully solving a vast majority of the instances (all in our case) even if at the expense of a large budget of function evaluations. Direct multisearch was applied off-the-shelf to determine the cardinality/mean-variance efficient frontier. The structure of problem [START_REF] Bertsimas | Algorithm for cardinality-constrained quadratic optimization[END_REF], or of its practical counterpart [START_REF] Cesarone | A new method for mean-variance portfolio optimization with cardinality constraints[END_REF], was essentially ignored. One can use the fact that the first objective function is smooth and of known derivatives to speed up the optimization and reduce even further the budget of function evaluations. Moreover, we also point out that it is trivial to run the poll step of direct multisearch in a parallel mode.

The use of derivative-free single or multiobjective optimization opens the research range of future work in sparse or dense portfolio selection. In fact, since derivative-free algorithms only rely on zero order information, they are applicable to any objective function of black-box type. One can thus use any measure to quantify the profit and risk of a portfolio. The classical Markowitz model assumes that the return of a portfolio is a linear combination of the returns of the individual securities. Also, it implicitly assumes a Gaussian distribution for the return, letting its variance be a natural measure of risk. However, it is known from the analysis of stylized facts that the distribution for the return of securities exhibits tails which are fatter than the Gaussian ones. Practitioners consider other measures of risk and profit better tailored to reality. Our approach to compute the cardinality/mean-variance efficient frontier is ready for application in such general scenarios.

A Using direct multisearch to determine efficient cardinality/mean-variance portfolios

A few modifications to problem [START_REF] Bertsimas | Algorithm for cardinality-constrained quadratic optimization[END_REF] were required to make it solvable by a multiobjective derivative-free solver, in particular by a direct multisearch one. In practice the first modification to (5) consisted of approximating the true cardinality, by introducing a tolerance , min

w∈R N -µ w w Qw min w∈R N N i=1 1 1 {|wi|> }
subject to e w = 1, N -1

L i ≤ w i ≤ U i , i =
i=1 1 1 {|wi|> } subject to -b ≤ w i ≤ b, i = 1, . . . , N -1, -b ≤ 1 - N -1 i=1 w i ≤ b, (10)
where w N in -µ w/w Qw was replaced by 1 -N -1 i=1 w i . We used all the default parameters of dms (version 0.2) with the following four exceptions. First, we needed to increase the maximum number of function evaluations allowed (from 20000 to 2000000 for N (= n) up to 50) given the dimension of our portfolios, as well as to require more accuracy by reducing the step size tolerance from 10 -3 to 10 -7 . Then we turned off the use of the cache of previously evaluated points to make the runs faster (the default version of dms keeps such a list to avoid evaluating points too close to those already evaluated). Lastly, we realized that initializing the list of feasible nondominated points with a singleton led to better results than initializing it with a set of roughly N points as it happens by default. Thus, we set the option list of dms to zero, which, given the bounds on the variables, assigns the origin to the initial list.

  min

  For the FTSE 100 data sets this sample is from 01/2003 to 12/2006 and for the FF data sets is from 07/1971 to 06/1996. Figures1-3 and Figures 4-6contain the plots of the efficient frontiers calculated for, respectively, the FTSE 100 and FF data sets. In all these plots we also marked three other portfolios. The first one is the 1/N portfolio corresponding to the naive strategy. A second one is obtained maximizing expected return per variance min to the extreme point (of maximum cardinality) of the efficient frontier (or Pareto front) of the cardinality/mean-variance biobjective optimization problem[START_REF] Bertsimas | Algorithm for cardinality-constrained quadratic optimization[END_REF]. The third one is a classical Markowitz related portfolio and is obtained by minimizing variance under no short-selling min w∈R N w Qw subject to e w = 1, w ≥ 0.

  By introducing binary variables, one can rewrite this problem as a mixed-integer quadratic programming (MIQP) problem: min w,y∈R N w Qw subject to e y ≤ K, e w = 1,

Figures 19 - 4 . 7

 1947 21 for the FTSE 100 portfolios and in Figures 22-24 for the FF ones. Results including the financial crisis years 2008-2010 The FTSE 100 data set used covered the period 2003-2007. With the aim of testing the robustness of the results, we also tried a FTSE 100 data set that covers the time window 2003-2010 (including thus the financial crisis years 2008-2010

1 ,subject to e w = 1 , 1 -

 111 . . . , N. chosen as = 10 -8 (1 1 represents the indicator function). Secondly, we selected symmetric bounds on the variables L i = -b and U i = b, -b ≤ w i ≤ b, i = 1, . . . , N, setting b = 10. Finally, we eliminated the constraint e w = 1 since direct search methods do not cope well with equality constraints. The version fed to the dms solver was then min w(1:N -1)∈R N -

Fig. 1 .Fig. 2 .

 12 Fig. 1. Efficient frontier of the biobjective cardinality/mean-variance problem for DTS1. Naive Markowitz mean per variance Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 3 .

 3 Fig. 3. Efficient frontier of the biobjective cardinality/mean-variance problem for DTS3. See the caption of Figure 1 for an explanation of the various symbols.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Efficient frontier of the biobjective cardinality/mean-variance problem for FF10. Naive Markowitz mean per variance Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 6 .

 6 Fig. 6. Efficient frontier of the biobjective cardinality/mean-variance problem for FF48. See the caption of Figure 4 for an explanation of the various symbols.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Out-of-sample performance for DTS1 measured by the Sharpe ratio over all the out-of-sample periods. --Naive -Markowitz mean per variance -.-Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 9 .

 9 Fig. 9. Out-of-sample performance for DTS3 measured by the Sharpe ratio over all the out-of-sample periods. See the caption of Figure 7 for an explanation of the various symbols and lines.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Out-of-sample performance for FF10 measured by the Sharpe ratio over all the out-of-sample periods. --Naive -Markowitz mean per variance -.-Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 12 .

 12 Fig. 12. Out-of-sample performance for FF48 measured by the Sharpe ratio over all the out-of-sample periods. See the caption of Figure 10 for an explanation of the various symbols and lines.

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Transaction costs of the efficient cardinality/mean-variance portfolios for DTS1. -Markowitz mean per variance -.-Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 15 .

 15 Fig. 15. Transaction costs of the efficient cardinality/mean-variance portfolios for DTS3. See the caption of Figure 13 for an explanation of the various symbols and lines.

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Transaction costs of the efficient cardinality/mean-variance portfolios for FF10. -Markowitz mean per variance -.-Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 18 .

 18 Fig. 18. Transaction costs of the efficient cardinality/mean-variance portfolios for FF48. See the caption of Figure 16 for an explanation of the various symbols and lines.

Fig. 19 .Fig. 20 .

 1920 Fig. 19. Out-of-sample performance for DTS1 measured by the Sharpe ratio of returns net of transaction costs over all the out-of-sample periods. --Naive -Markowitz mean per variance -.-Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 21 .

 21 Fig. 21. Out-of-sample performance for DTS3 measured by the Sharpe ratio of returns net of transaction costs over all the out-of-sample periods. See the caption of Figure 19 for an explanation of the various symbols and lines.

Fig. 22 .

 22 Fig. 22. Out-of-sample performance for FF10 measured by the Sharpe ratio of returns net of transaction costs over all the out-of-sample periods. --Naive -Markowitz mean per variance -.-Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 24 .

 24 Fig. 24. Out-of-sample performance for FF48 measured by the Sharpe ratio of returns net of transaction costs over all the out-of-sample periods. See the caption of Figure 22 for an explanation of the various symbols and lines.

Fig. 25 .Fig. 26 .

 2526 Fig. 25. Out-of-sample performance for DTS1, including the financial crisis years 2008-2010, measured by the Sharpe ratio over all the out-of-sample periods. --Naive -Markowitz mean per variance -.-Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 27 .

 27 Fig. 27. Out-of-sample performance for DTS3, including the financial crisis years 2008-2010, measured by the Sharpe ratio over all the out-of-sample periods. See the caption of Figure 25 for an explanation of the various symbols and lines.

Fig. 28 .

 28 Fig. 28. Transaction costs of the efficient cardinality/mean-variance portfolios for DTS1, including the financial crisis years 2008-2010. -Markowitz mean per variance -.-Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 30 .

 30 Fig. 30. Transaction costs of the efficient cardinality/mean-variance portfolios for DTS3, including the financial crisis years 2008-2010. See the caption of Figure 28 for an explanation of the various symbols and lines.

Fig. 31 .

 31 Fig. 31. Out-of-sample performance for DTS1, including the financial crisis years 2008-2010, measured by the Sharpe ratio of returns net of transaction costs over all the out-of-sample periods. --Naive -Markowitz mean per variance -.-Markowitz minimum variance • cardinality/mean-variance cardinality constrained minimum variance

Fig. 33 .

 33 Fig. 33. Out-of-sample performance for DTS3, including the financial crisis years 2008-2010, measured by the Sharpe ratio of returns net of transaction costs over all the out-of-sample periods. See the caption of Figure 31 for an explanation of the various symbols and lines.

Fig. 34 .

 34 Fig. 34. Efficient frontier of the biobjective cardinality/mean-variance problem for FF100.

For the FTSE 100 data sets we used the 3 month Treasury-Bills UK. Such data is public and made available by the Bank of England, at the site http://www. bankofengland.co.uk. For the FF data sets we used the 90-day Treasury-Bills US. Such data is public and made available by the Federal Reserve, at the site http://www.federalreserve.gov.

We used as a risk-free asset the daily startling certificate of deposit interest rate. Such data is public and made available by the Bank of England, at the site http: //www.bankofengland.co.uk.

Support for this research was provided by FCT under the grant PTDC/MAT/098214/2008.


Simone Cacace

Emiliano Cristiani

Maurizio Falcone

Two Semi-Lagrangian Fast Methods for Hamilton-Jacobi-Bellman Equations

Keywords: Single-pass methods, fast iterative method, fast sweeping method, fast marching method

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The study of Hamilton-Jacobi (HJ) equations arises in several contexts, including classical mechanics, front propagation, control problems and differential games. In particular, for optimal control problems, the value function can be characterized as the unique viscosity solution of a Hamilton-Jacobi-Bellman (HJB) equation. Unfortunately, solving numerically the HJB equation can be rather expensive from the computational point of view. This is the reason why in the last years an increasing number of efficient techniques have been proposed, see, e.g., [START_REF] Cacace | Can local single-pass methods solve any stationary Hamilton-Jacobi-Bellman equation[END_REF] for a brief review.

Basically, these algorithms are divided in two main classes: single-pass and iterative. An algorithm is said to be single-pass if one can fix a priori a (small) number r which depends only on the equation and on the mesh structure (not on the number of mesh points) such that each mesh point is re-computed at most r times. Single-pass algorithms usually divide the numerical grid in, at least, three changing subsets: Accepted (ACC), Considered (CONS), and Far (FAR). Nodes in ACC are definitively computed, nodes in CONS are computed but their values are not yet final, and nodes in FAR are not yet computed. We say that a single-pass algorithm is local if the computation at any mesh point involves only the values of first neighboring nodes, the region CONS is 1-cell-thick and no information coming from FAR region is used. The methods which are not single-pass are iterative.

Among fast methods, the prototype algorithm for the local single-pass class is the Fast Marching Method (FMM) [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF], while that for the iterative class is the Fast Sweeping Method (FSM) [START_REF] Kao | Lax-Friedrichs sweeping scheme for static Hamilton-Jacobi equations[END_REF][START_REF] Qian | A fast sweeping method for static convex Hamilton-Jacobi equations[END_REF][START_REF] Tsai | Fast sweeping algorithms for a class of Hamilton-Jacobi equations[END_REF][START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF]. Another interesting method is the Fast Iterative Method (FIM) [START_REF] Fu | A fast iterative method for solving the eikonal equation on triangulated surfaces[END_REF][START_REF] Jeong | A fast iterative method for eikonal equations[END_REF][START_REF] Jeong | A fast iterative method for a class of Hamilton-Jacobi equations on parallel systems[END_REF], which shares some features with both iterative and single-pass methods. Recently, Cacace et al. [START_REF] Cacace | Can local single-pass methods solve any stationary Hamilton-Jacobi-Bellman equation[END_REF] have shown that it is not possible to create a local single-pass algorithm for solving general HJB equations. This motivates the efforts to develop new techniques, particularly in the class of iterative methods.

In this paper, we focus on the following minimum time HJB equation

sup a∈B1 {-f (x, a) • ∇T (x)} = 1, x ∈ R d \T ( 1
)
where To simplify the notations, we restrict the discussion to the case d = 2. Generalizations of the considered algorithms to any space dimension is straightforward, although the implementation is not trivial. The goal of this paper is twofold: First, we investigate the possibility of applying a semi-Lagrangian version of the FIM to equation (1). To our knowledge, FIM was only used for solving the eikonal equation [START_REF] Jeong | A fast iterative method for eikonal equations[END_REF] and a special class of HJ equations [START_REF] Jeong | A fast iterative method for a class of Hamilton-Jacobi equations on parallel systems[END_REF], although there is no particular constraint to apply it in a more general framework. The algorithm indeed does not rely on the special form and features of the eikonal equation. In addition, we measure the degree of "iterativeness" of FIM, keeping track of how many times each grid node is inserted into the list of nodes which are actually computed at each step. Interestingly, the results indirectly confirm the findings of [START_REF] Cacace | Can local single-pass methods solve any stationary Hamilton-Jacobi-Bellman equation[END_REF], showing that general HJB equations require the nodes to be visited an (a priori ) unknown number of times, i.e. single-pass methods do not apply.

d is the space dimension, T is a closed nonempty target set in R d , f : R d × B 1 → R d is
Second, we propose a new acceleration technique for the FSM, which is effective when (1) is discretized by means of a semi-Lagrangian scheme (see [START_REF] Falcone | Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations[END_REF] for a comprehensive introduction). It reduces the CPU load for the sup search in (1), neglecting the control directions which are downwind with respect to the current sweep. The new method results to be remarkably faster, although (in general) the number of iterations needed for convergence increases.

Semi-Lagrangian approximation

Let us introduce a structured grid G and denote its nodes by x i , i = 1, . . . , N . The space step is assumed to be uniform and equal to ∆x > 0. Standard argu-ments [START_REF] Falcone | Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations[END_REF] lead to the following discrete version of equation (1):

T (x i ) ≈ T (x i ) = min a∈B1 T (x i,a ) + |x i -xi,a | |f (x i , a)| , x i ∈ G (2)
where xi,a is a non-mesh point, obtained by integrating, until a certain final time τ , the ordinary differential equation

ẏ(t) = f (y, a), t ∈ [0, τ ] y(0) = x i (3)
and then setting xi,a = y(τ ). To make the scheme fully discrete, the set of admissible controls B 1 is discretized with N c points and we denote by a * the optimal control achieving the minimum. Note that we can get different versions of the semi-Lagrangian (SL) scheme (2) varying τ , the method used to solve (3), and the interpolation method used to compute T (x i,a ). Moreover, we remark that, in any single-pass method, the computation of T (x i ) cannot involve the value T (x i ) itself, because this self-dependency would make the method iterative.

Here we use a 3-point scheme: Equation ( 3) is solved by an explicit forward Euler scheme until the solution is at distance ∆x from x i , where it falls inside the triangle of vertices x i,1 , x i,2 , and x i,3 , to be chosen among the first neighbors of x i . The value T (x i ) is computed by a two-dimensional linear interpolation of the values T (x i,1 ), T (x i,2 ) and T (x i,3 ) (see [START_REF] Cacace | Can local single-pass methods solve any stationary Hamilton-Jacobi-Bellman equation[END_REF] for details).

Limits of local single-pass methods

In this section we briefly recall the main result of [START_REF] Cacace | Can local single-pass methods solve any stationary Hamilton-Jacobi-Bellman equation[END_REF]. From the numerical point of view, it is meaningful to divide HJB equations into four classes. For any given mesh, we have:

(ISO) Equations whose characteristic curves coincide or lie in the same simplex of the gradient curves of their solutions. (¬ISO) Equations for which there exists at least a grid node where the characteristic curve and the gradient curve of the solution do not lie in the same simplex.

(REG) Equations with non-crossing (regular) characteristic curves. Characteristics spread from the target T to the rest of the domain without intersecting. (¬REG) Equations with crossing characteristic curves. Characteristics start from the target T and then meet in finite time, creating shocks. As a result, the solution T is not differentiable at shocks.

Let us summarize here the main remarks on single-pass methods: • FMM works for equations of type ISO and fails for equations of type ¬ISO (see [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: theory and algorithms[END_REF] for further details and explanations), while FSM can be successfully applied in any case.

• Handling ¬ISO case requires CONS not to follow the level sets of the solution itself. Indeed, if CONS turns out to be an approximation of the level sets of the solution, it means that the solution is computed in an increasing order, thus following the gradient curve rather than the characteristic curve.

• Handling ¬REG case requires CONS to be an approximation of the level sets of the solution. Let us clarify this point. Consider the ¬REG case and let x be a point belonging to a shock, i.e. where the solution is not differentiable. By definition, the value T (x) is carried by two or more characteristic curves reaching x at the same time. Similarly, let x i be a grid node ∆x-close to the shock. In order to mimic the continuous case, x i has to be approached by the ACC region approximately at the same time from the directions corresponding to the characteristic curves. In this case, the value T (x i ) is correct (no matter which upwind direction is chosen) and, more important, the characteristic information stops at x i and it is no longer propagated, getting stuck by the ACC region. As a consequence, the shock is localized properly. On the other hand, if CONS region is not an approximation of a level set of the solution a node x i close to a shock can be reached by ACC at different times. When ACC reaches x i for the first time, it is impossible to detect the presence of the shock by using only local information. Indeed, only a global view of the solution allows one to know that another characteristic curve will reach x i at a later time. As a consequence, the algorithm continues the enlargement of CONS and ACC, thus making an error that cannot be corrected by the following iterations.

In conclusion, we get that local single-pass methods cannot handle equations ¬ISO & ¬REG. In this situation, one has to add non local information regarding the location of the shock, or going back to nodes in ACC at later time, breaking the single-pass property. This motivates the investigation of new techniques, especially iterative methods, as the ones described in the next sections.

Fast iterative method

In this section we briefly recall the construction of FIM [START_REF] Fu | A fast iterative method for solving the eikonal equation on triangulated surfaces[END_REF][START_REF] Jeong | A fast iterative method for eikonal equations[END_REF][START_REF] Jeong | A fast iterative method for a class of Hamilton-Jacobi equations on parallel systems[END_REF]. As in FMM, the main idea of FIM is to update only few grid nodes at each step. These nodes are stored in a separated list, called active list. During each step, the list of active nodes is modified, and the band thickens or expands to include all nodes that could be affected by the current updates. A node can be removed from the active list when its value is up to date with respect to its neighbors (i.e., it has reached convergence) and can be appended to the list (relisted) whenever any upwind neighbor's value has changed.

FIM is formally an iterative method, since the number of times a grid node is visited depends on the dynamics and on the grid size. On the other hand, the active list resembles the set CONS of FMM, and in some special cases FIM is in fact a single-pass algorithm, see Section 6. Nevertheless, the active list and CONS differ for some important features. The first is that the active list is not kept ordered, and then the causality relationship among grid nodes is lost. The second is that the active list can be more than 1-node-thick, i.e. it can approximate a two-dimensional set. Finally, grid nodes removed from the active list can re-enter at a later time. This is the price to pay for loosing the causality.

The FIM algorithm consists of two parts, the initialization and the updating. In the initialization step, one has to set the boundary conditions and set the values of the rest of the grid nodes to infinity (or some very large value). Next, the adjacent neighbors of the source nodes (i.e. the target) are added to the active list. In the updating step, for every point in the list, one computes the new value and checks if the value at the node has converged by comparing the old and the new value at the considered point. If it has converged, one removes the node from the list and append to the list any non active adjacent node such that its updated value is less than the current one. The algorithm runs until the list is empty.

FIM was introduced for solving a special class of HJ equations [START_REF] Jeong | A fast iterative method for eikonal equations[END_REF][START_REF] Jeong | A fast iterative method for a class of Hamilton-Jacobi equations on parallel systems[END_REF]. Nevertheless, in Section 6 we show that FIM based on a SL discretization can be successfully applied to general HJB equations with no modifications.

An optimized fast sweeping method

FSM is another popular method for solving HJ equations [START_REF] Kao | Lax-Friedrichs sweeping scheme for static Hamilton-Jacobi equations[END_REF][START_REF] Qian | A fast sweeping method for static convex Hamilton-Jacobi equations[END_REF][START_REF] Tsai | Fast sweeping algorithms for a class of Hamilton-Jacobi equations[END_REF][START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF]. The main advantage of the method is its implementation, which is extremely easy (easier than that of FMM and FIM). FSM is basically the classical iterative (fixed-point) method, since each node is visited in a predefined order, until convergence is reached. Here, the visiting directions (sweeps) are alternated in order to follow all possible characteristic directions, trying to exploit causality. In two-dimensional problems, the grid is visited sweeping in four directions:

S → N & W → E, S → N & E → W , N → S & E → W and N → S & W → E.
The key point is the Gauss-Seidel-like update of grid nodes, which allows one to compute in a cascade fashion a relevant part of the grid nodes in only one sweep. Indeed it is well known that in the case of eikonal equations FSM converges in only four sweeps [START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF].

Here we propose an easy modification of the FSM based on a SL discretization, aiming at saving CPU time for each sweep. Let us explain the idea in the case of a dynamics of the form c(x, a)a, with c > 0. It is clear that during the sweep S → N & W → E the algorithm cannot exploit the power of the Gauss-Seidel cascade for the information coming from NE. Indeed, even if a node actually depends on its NE neighbor, that information flows upwind and it is not propagated to other nodes during the current sweep. Then, we propose to remove downwind discrete controls from the minimum search in the SL scheme (2), since they have small or no effect in the update of the nodes, see Figure 1.

The assumption c > 0 is needed to preserve the order of the quadrants between the control a and the resulting dynamics c(x, a)a. Otherwise, the choice of controls to be removed should be adapted according to the sign of c. Note that the control set B 1 is reduced to a upwind 3/4 of ball. Then, let us denote by Upwind Fast Sweeping Method 3/4 (UFSM3/4) the classical FSM with this control reduction. An additional speedup, that we expect to work only in case where the characteristics are essentially straight, consists in reducing further the control ball to a upwind 1/4 of ball (UFSM1/4).

Numerical experiments

In this section we compare the performance of FSM, FIM, UFSM1/4 and UFSM3/4 on the following equations:

Equation Dynamics Class HJB-1 f (x, y, a) = a ISO & REG HJB-2 f (x, y, a) = (1 + 4χ {x>1} ) a ISO & ¬REG HJB-3 f (x, y, a) = m λ,µ (a) a ¬ISO & REG HJB-4 f (x, y, a) = F 2 (x, y)m p(x,y),q(x,y) (a) a ¬ISO & ¬REG HJB-5 f (x, y, a) = (1 + |x + y|)m λ,µ (a) a ¬ISO & ¬REG
where we defined m λ,µ (a) = (1 + (λ a 1 + µ a 2 ) 2 ) -1 2 for λ, µ ∈ R and we denoted by χ S the characteristic function of a set S. Moreover, for c 1 , c 2 , c 3 , c 4 > 0, we defined

C(x) = c 1 sin c 2 πx c 3 + c 4 , F 1 (x, y), F 2 (x, y) = (0.5, 1) if y ≤ C(x) (2, 3) otherwise , M (x, y) = F 2 2 (x,y) F 2 1 (x,y) -1 1 + C 2 (x)
, p(x, y) = M (x, y)C (x), q(x, y) = -M (x, y) .

In all the following tests we set Ω = [-2, 2] 2 (except Test 4), the target T = (0, 0) and the number of discrete controls N c = 32. Regarding FIM, we keep track of the history of the active list by counting the number I i of times the node x i enters the active list. The number I max := max i I i gives a measure of the "iterativeness" of the method.

Test 1. Here we solve equation HJB-1. Figure 2 shows the evolution of FIM's active list. Differently from FMM, where the CONS set expands from the target following concentric circles (i.e. the level sets of the solution), here the active set moves following concentric squares (cf. the behavior of the CONS region of the safe method studied in [START_REF] Cacace | Can local single-pass methods solve any stationary Hamilton-Jacobi-Bellman equation[END_REF]). As one can expect I max = 1, meaning that FIM behaves like a single-pass method. Table 1 compares CPU times of the methods on different grids and the number of sweeps needed by sweeping methods to reach convergence. FSM converges in 4 sweeps for this equation, the additional sweep Fig. 2. FIM for HJB-1: active nodes at different steps reported in Table 1 is the one required by the algorithm to check convergence. All the methods compute the same solution. In particular we see that FIM is slightly slower than FSM, as noted in [START_REF] Jeong | A fast iterative method for eikonal equations[END_REF]. On the other hand, UFSM methods (both 1/4 and 3/4) still converge in 5 sweeps, thus overcoming FSM.

Test 2. Here we solve equation HJB-2. Figure 3 shows the optimal vector field f (x, a * ) and the history of active nodes (in grey scale, where black corresponds to I max and white to 0). The maximal number of re-activation is I max = 3 and re-activation of nodes appears for the first time close to the shock line, see Figure 3-center. This depends on the fact that the active list is not an approximation of a level set of the solution and the equation HJB-2 falls in the ¬REG class. Then FIM is not able to capture the shock properly (see Section 3) in a single-pass fashion, but has to come back to recompute wrong values. Table 1 compares the methods. Results are similar to those of the previous test. Test 3. Here we solve equation HJB-3 for λ = 10 and µ = 5, namely the anisotropic eikonal equation, a well known example where FMM fails in computing the correct solution, due to the fact that characteristics do not coincide with gradient curves of the solution, see [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: theory and algorithms[END_REF]. In this case FIM let evolve the active list as for the eikonal equation (Test 1, Figure 2) and produces a maximal number of re-activation I max = 1. Again, this means that equation HJB-3 can be successfully solved by a local single-pass method, as the safe method introduced in [START_REF] Cacace | Can local single-pass methods solve any stationary Hamilton-Jacobi-Bellman equation[END_REF] for the class REG. We refer to Table 1 for a comparison of the methods. All the methods compute the same solution, meaning that FIM can work for equations with substantial anisotropy and inhomogeneities (see also next test). Unexpectedly, also UFSM1/4 is able to correctly follow quite curved characteristics, see Figure 4-left). Results in Table 1 show that sweeping methods need a large number of sweeps to reach convergence (even more for UFSMs, due to the control set reduction). This makes FIM be the fastest method. The maximal number of re-activation for the active list is I max = 7 and Figure 4-center/right shows that re-activation of nodes appears both close to the shocks and where the optimal field exhibits rapid changes of direction. Test 5. Here we solve HJB-5 for λ = 10 and µ = 5. This is the hardest example of class ¬ISO & ¬REG presented in [START_REF] Cacace | Can local single-pass methods solve any stationary Hamilton-Jacobi-Bellman equation[END_REF], where the shock (see the cubic-like curve in Figure 5-left/center) and a strong anisotropy region meet at the target. Sweeping methods FSM and UFSM3/4 require much more sweeps with respect to the previous tests, while UFSM1/4 fails in computing the correct solution, confirming that the control set reduction to 1/4 of ball cannot be applied in any case.

The maximal number of re-activation for FIM is I max = 30 (see Figure 5right), whereas the evolution of the active list is extremely complicated and also produces regions of dimension two (see Figure 6). Nevertheless, results in Table 1 shows that, as the grid increases, FIM is still the fastest method. The presence of a two-dimensional active list clearly proves that local single-pass methods cannot be applied since the enlargement of CONS is required (cf. the buffered fast marching method [START_REF] Cristiani | A fast marching method for Hamilton-Jacobi equations modeling monotone front propagations[END_REF]).

Conclusions

Tests performed in Section 6 show that FIM can be successfully used for solving general HJB equations with no modifications with respect to the original algorithm. Moreover, FIM appears to be the fastest method in case of complicated ¬ISO & ¬REG equations because of the large number of iterations needed by the sweeping methods. Considering that the implementation of FIM is not harder than that of FSM, we think that, overall, FIM is the best method among the tested ones.

UFSM3/4 is always preferable to FSM, since the larger number of iterations needed for convergence are widely counterbalanced by the speedup for each single sweep. UFSM1/4 is instead not safely applicable for general equations.

Finally, the results of this paper confirm those of [START_REF] Cacace | Can local single-pass methods solve any stationary Hamilton-Jacobi-Bellman equation[END_REF]. Complicated ¬ISO & ¬REG equations require to pass through some nodes more than one time (cf., e.g., Figure 5-right and [1, Figure 7]), and exhibit two-dimensional regions in which every node depends on each other.

  a given vector-valued Lipschitz continuous function, and B 1 is the unit ball in R d centered in the origin, representing the set of the admissible controls. We complement the equation with homogeneous Dirichlet condition T = 0 on T . Let us note that if f (x, a) = c(x)a, equation (1) becomes the eikonal equation c(x)|∇T (x)| = 1.

Fig. 1 .

 1 Fig.1. Downwind controls with respect to the four sweeps: dashed arcs identify the directions to be removed

0 3 Fig. 3 .

 33 Fig. 3. FIM for HJB-2: Optimal controls (left), re-activation history at an intermediate step (center), re-activation history and level sets of the solution (right)

Test 4 .

 4 Here we solve equation HJB-4 in Ω = [-0.5, 0.5] 2 for c 1 = 0.1225, c 2 = 2, c 3 = 0.5 and c 4 = 0, an example of class ¬ISO & ¬REG coming from seismic imaging. It is a inhomogeneous anisotropic eikonal equation on a domain with two layers separated by a sinusoidal profile C(x), with different constant anisotropy coefficients in each layer (given by the pairs (F 1 , F 2 ) = (0.5, 1) and (F 1 , F 2 ) = (2, 3)).

0 7 Fig. 4 .

 74 Fig. 4. FIM for HJB-4: Optimal controls (left), re-activation history at an intermediate step (center), re-activation history and level sets of the solution (right)

Fig. 5 .Fig. 6 .

 56 Fig. 5. FIM for HJB-5: Optimal controls (left), level sets of the solution (center), reactivation history (right)

Table 1 .

 1 CPU times (seconds) and number of sweeps. Fastest method in bold

	Equation Grid ∆x FSM (sweeps) FIM UFSM1/4 (sweeps) UFSM3/4 (sweeps)
	HJB-1 101 0.04	0.13 (5)	0.17	0.04 (5)	0.10 (5)
	HJB-1 201 0.02	0.51 (5)	0.72	0.15 (5)	0.40 (5)
	HJB-1 401 0.01	2.04 (5)	3.21	0.58 (5)	1.51 (5)
	HJB-2 101 0.04	0.16 (5)	0.21	0.04 (5)	0.12 (5)
	HJB-2 201 0.02	0.63 (5)	0.87	0.19 (5)	0.46 (5)
	HJB-2 401 0.01	2.46 (5)	3.80	0.70 (5)	1.84 (5)
	HJB-3 101 0.04	0.31 (5)	0.38	0.09 (5)	0.23 (5)
	HJB-3 201 0.02	1.23 (5)	1.56	0.35 (5)	0.88 (5)
	HJB-3 401 0.01	4.88 (5)	6.57	1.38 (5)	3.53 (5)
	HJB-4 101 0.01	5.72 (25)	1.93	2.18 (34)	5.62 (34)
	HJB-4 201 0.005 22.70 (25)	7.68	8.66 (34)	19.58 (30)
	HJB-4 401 0.0025 99.38 (28) 29.36	34.07 (34)	77.14 (30)
	HJB-5 101 0.04	3.23 (53)	5.30	-	2.38 (53)
	HJB-5 201 0.02	13.30 (55)	3.32	-	9.78 (55)
	HJB-5 401 0.01	52.93 (55) 14.53	-	39.16 (55)

This research was supported by the following grants: AFOSR Grant FA9550-10-1-0029, ITN-Marie Curie Grant 264735-SADCO.