Guillaume Biessy

Smoothing of multidimensional biometric laws in a Long-Term Care Insurance portfolio

, which brings an additional dimension to the estimation problem. In French private LTC insurance, products are purchased at the age of 60 while most claims occur after 85, more than 25 years later, making it an actual Long-Term risk alongside Longevity. Therefore, it is paramount for insurers to develop a deep knowledge of the biometric laws in order to get a clear vision of their liabilities and increase insurance premiums if necessary.

 to smooth the obtained rates and derive extrapolations. Our modelling approach naturally fits in the framework introduced in Currie (2013) for inference of stochastic models with P-splines and constraints. Our implementation further capitalize on the multidimensional structure of our data through the use of Generalized Linear Array Models to greatly improve computation times. We find the restricted interaction approach to be both accurate and parsimonious with the additional advantage of having a clearer interpretation of the results. Our approach could easily be applied to other insurance risks such as disability or critical illness.

Introduction

Over the past decades, due to progress in medicine and especially in the treatment of cardiovascular diseases, there has been a sustained increase in life expectancy, especially at higher ages, leading to a tremendous increase in the number of elderly people. Among other consequences, this has considerably increased the need for Long-Term Care (LTC) services. As LTC is very expansive when compared to the average earning of retired elderly people, this led to the launch of the first LTC products in the late 1970s in the US and in 1985 for France.

While LTC has received a lot of attention from both economists and the insurance industry, publications on LTC modelling remain very scarce, despite strong similarities with epidemiology where multi-state processes are commonly used as well as demography for problems related to old age mortality extrapolation and forecasting.

Among the 3 biometric laws, mortality in LTC proves difficult to model because LTC can be caused by various pathologies such as cancer, dementia, neurological or cardiovascular diseases. All those pathologies are linked to ageing but they are very different, resulting in high heterogeneity in trajectories of disabled people. This leads to a mixture effect that has been identified as the cause for a very high mortality during the first months following entry in the LTC state [START_REF] Biessy | Semi-Markov modeling of the loss of autonomy among elderly people: application to long-term care insurance[END_REF]. Semi-Markov models where the time spent in the LTC state is used as a predictor for the mortality, have been identified for long as the adequate tool to handle this phenomenon.

Nevertheless attempts to derive 2-dimension estimates of mortality in LTC remain scarce. [START_REF] Tomas | Multidimensional smoothing by adaptive local kernel-weighted log-likelihood: Application to long-term care insurance[END_REF] provides a methodology based on local-likelihood smoothing, working on monthly mortality rates for which Chapter 4 of [START_REF] Biessy | Semi-Markov modeling of the loss of autonomy among elderly people: application to long-term care insurance[END_REF] provides a continuous-time variant. [START_REF] Biessy | Continuous-time semi-markov inference of biometric laws associated with a long-term care insurance portfolio[END_REF] introduces a parametric model to replicate the mixture effect observed in the data. [START_REF] Guibert | Non-parametric inference of transition probabilities based on aalen-johansen integral estimators for acyclic multi-state models: application to ltc insurance[END_REF] introduce non-parametric estimators based on Aalen-Johansen integrals.

On the other hand, due to the growing stress put on pension funds by the increase in life expectancy, sophisticated multidimensional models have emerged to estimate and forecast mortality in order to support the booming business of Longevity transaction between reinsurers and pension funds. Among the first prospective models are the Age-Period, Age-Cohort and Age-Period-Cohort model introduced in Clayton and Schifflers (1987a) and Clayton and Schifflers (1987b). Then follows the widespread model of [START_REF] Lee | Modeling and forecasting US mortality[END_REF] further improved when [START_REF] Brouhns | A poisson log-bilinear regression approach to the construction of projected lifetables[END_REF] suggested to work on the number of deaths rather than the mortality rates which led to a more convenient formulation of the original model. [START_REF] Currie | Smoothing and forecasting mortality rates[END_REF] suggested an alternative approach based on multidimensional smoothing using P-splines [START_REF] Eilers | Flexible smoothing with B-splines and penalties[END_REF] for which the Generalized Linear Array Model (GLAM) framework introduced in [START_REF] Currie | Generalized linear array models with applications to multidimensional smoothing[END_REF] provided much faster computation. Finally [START_REF] Currie | Smoothing constrained generalized linear models with an application to the lee-carter model[END_REF] gave a more general framework for the fitting of the Lee-Carter model which could be extended to fit even more sophisticated stochastic models with cohort effects.

To model Long-Term Care we rely on the illness-death model which was first introduced in [START_REF] Kodell | An illness-death model for the study of the carcinogenic process using survival/sacrifice data[END_REF]. We make the additional assumption that no transition from the LTC state back to the Healthy state may occur. Indeed, due to the definition of LTC in France, such recoveries are very unlikely to occur. Our data support this assumption with fewer than one case of recovery being reported per thousands of claims. Finally, we focus on biometric risks and do not model the lapse risk. Figure 1 give a summary of the model where we added the transition toward lapsed for completeness. We therefore model 3 transitions: active lives may move to LTC or death so we have two competing risks and disabled lives may only move to the death state. We allow all transitions to depend on age x and gender g. We let z denote the time spent in the current state and assume that the mortality in LTC may depend on z where time spent in autonomy is assumed to have no effect on either of the 3 transitions. The final model may be seen as a particular case of a semi-Markov model.

Healthy LTC

Lapsed Death

λ 01 (x, g) λ 02 (x, g) λ 12 (x, z, g)
Figure 1: The Long-Term Care model with one level of care and no recoveries from LTC

In this paper, we present a new methodology, inspired by the aforementioned papers, to infer those 3 biometric laws. We rely on P-spline for smooth inference of each biometric law as well as extrapolation for higher ages. In the case of mortality in LTC, we account for time spent in LTC and develop two different approaches to smooth the mortality surface.

The first section illustrates the P-spline methodology that is key to our inference procedure. The second section introduces the Poisson GLM framework in which our models for each of the biometric laws are defined, as well as the scoring algorithm that will be used to fit the parameters. Finally it provides some computational details of our approach, especially the GLAM framework presented in [START_REF] Currie | Generalized linear array models with applications to multidimensional smoothing[END_REF] which allows for faster computation in the case of multidimensional problems. The third section provides smooth estimates of autonomous mortality and incidence in LTC with their extrapolations at higher ages. The penultimate section consider the case of mortality in LTC for which two different approaches are developed and then compared: unconstrained smoothing in 2 dimensions on one hand and a decomposition over the age and duration in LTC component with restricted interactions on the other hand. Finally, we discuss the results and potential extensions to the proposed methodology.

P-splines smoothing

Our methodology relies heavily on P-splines smoothing which is incorporated in our inference procedure. The literature about P-splines is abundant and we refer the interested reader to [START_REF] Eilers | Flexible smoothing with B-splines and penalties[END_REF] for the founding paper. In what follows we merely highlight the main points in the methodology. P-splines stands for Penalized basis-splines. It combines basis splines (abbreviated in B-Splines) formally introduced in De Boor (1978) with penalized likelihood. Psplines is easy to both understand and implement, making it a very good all-around technique. Basis splines are piecewise polynomial functions of a given degree d that are non-zero only over d + 1 consecutive intervals, delimited by d knots. The B-splines as well as their first d -1 derivatives are continuous in each knot. With all those properties, it can be shown that the Bsplines are uniquely defined up to a scaling factor. Figure 2 shows a representation of B-splines of degrees 0 o 3. The interested reader may refer to De Boor (1978) for a very complete book about B-splines that also illustrates how they may be constructed recursively. In what follows, we consider equally spaced knots and splines of degree d = 3.

B-Splines

In order to perform smoothing using B-splines, the first step is to build a full B-splines basis. First the support for the observation must be divided into n d subdivisions delimited by n d + 1 equally-spaced knots called internal knots. To complete the basis, d external knots must also be added on each side, increasing the total number of knots to n d + 2d + 1. The B-splines basis then consists in n d + d splines with identical scaling factor. For each point within the range delimited by the internal knots, d + 1 splines are non-zero (only d at the positions of internal knots) and their sum is equal to 1. Figure 3 represents an example of such basis. Using the basis previously defined, we introduce the dummy function

f (x) = 2 + cos x -2 2
and assume our observations are sampled values of f (x) for x ∈ {3. 05, 3.15, . . . , 6.85, 6.95} with a white noise of standard deviation σ = 0.2. If θ represents the vector of weights for the B-splines, then the regression matrix B for the basis is such that Bθ is the vector of fitted values at the observation points. In the Gaussian case where y = Bθ + , ∼ N (0, σ 2 I n ) the fitted values ŷ are obtained from the observations y using θ = (B B) -1 B y and ŷ = B θ. An illustration of B-splines smoothing with the dummy function and the basis is provided by Figure 4.

Penalized likelihood

The idea of introducing a penalization term to smooth the data can be traced back to [START_REF] Whittaker | On a new method of graduation[END_REF]. If (θ) is the log-likelihood function for the model, the penalized log-likelihood function P (θ) is given by

P (θ) = (θ) + 1 2 θ P θ
where P is a penalty matrix.

For penalties over a single dimension, the penalty matrix P may be defined as P π = πD k,p D k,p where D k,p is a matrix of dimensions (p -k) × p is obtained by differentiating k times I p the identity matrix of size p where k ∈ N * , p is the length of θ and π is a scaling parameter. The matrices D 1,p and D 2,p are as follows:

D 1,p =       1 -1 0 . . . 0 0 0 1 -1 . . . 0 0 . . . . . . . . . . . . . . . 0 0 0 . . . 1 -1       D 2,p =       1 -2 1 0 . . . 0 0 0 0 1 -2 1 . . . 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 . . . 1 -2 1      
Refined penalty matrices may also be obtained by combining penalties matrix of different orders, for example:

P π = πD 2,p D 2,p + 2 √ πD 1,p D 1,p

P-Splines

In P-splines smoothing, B-splines and penalized likelihood are combined. In B-Splines smoothing, once the degree of the basis has been picked, the number of subdivisions is the only adjustment parameter. Increasing the number of splines results in a more accurate fit to the data with the risk of over-fitting. An optimal number of subdivisions may be determined using cross-validation. However, there may be areas in which a higher density of splines is required when other areas are very regular. A solution can be to use non-equally spaced knots but this comes to the detriment of the simplicity of the method. P-splines offer an attractive solution with the penalization term which takes care of avoiding over-fitting, no matter the number of splines. Therefore, this number of splines should only by picked on the basis of computational limitations. In P-splines smoothing, the order and scale of the penalty term becomes the main adjustment parameters. Higher order penalties result in smoother curves. However, the order of the penalty also plays an important role when it comes to extrapolation. Indeed, a penalty of order k will not penalize polynomial functions of up to degree k -1, meaning the asymptote for the extrapolation will be of degree k -1. A penalty order k = 2 may therefore be seen as a suitable default choice with a reasonable linear extrapolation while k = 1 may be conveniently used to force an horizontal asymptote to the fitted curve. The scale of the penalty π may be determined by minimizing a criterion such as the Bayesian Information Criterion (BIC): BIC = dev + log n × df where dev = 2 [ max -(θ)] is the deviance of the model (l being the log-likelihood function and l max the maximum reachable log-likelihood) and df = tr(H) the effective degrees of freedom (where n is the number of observations and H is the hat matrix for the estimation problem).

3 Estimation framework

Poisson GLM framework

The framework of Generalized Linear Models (GLM) was introduced by [START_REF] Nelder | Generalized linear models[END_REF] and unify the properties of various distribution of the exponential family. A GLM is defined by specifying 3 elements: a random component Y with mean µ = E(Y ), a deterministic component η = Xθ + offset and a link function g linking both previous components so that η = g(µ). For each family any bijective function may be chosen as link yet some link, called canonical links, functions are preferred because they make computations easier. In what follows we work on the random variable Y = D, the number of observed deaths (or entries in LTC when working on the transition from autonomous to LTC) and assume it follows a Poisson distribution of parameter λe c where λ is the transition intensity and e c the central exposure, the sum of all observed lifetimes. We use the canonical link for the Poisson distribution which is the log link.

Our model specification is hence given by log E(D) = log µ = η = log λ + log e c . We further note log λ = Xθ + ξ where ξ is the part of log λ that does not depend on the parameter θ then we can write log E(D) = Xθ + offset where offset = log e c + ξ. In what follows, every model we introduce may be fitted by solving one or two equations of this form.

The Poisson distribution assumption comes with a number of more and less obvious underlying assumptions. First to assume the decomposition µ = λe c , the transition intensity must be considered as constant over the interval where the central exposure is computed. We consider the transition intensity to be constant over one-year intervals of the form [x, x+1) for x ∈ N. For the study of mortality in LTC this assumption is too restrictive due to the high initial mortality and we consider instead one-month intervals for the first year in LTC. The Poisson distribution can also give an unbounded number of deaths while the observed number of deaths is actually limited by the number of individuals under risk. Then, by assuming a Poisson distribution of parameter λe c we implicitly consider the central exposure a deterministic quantity while it depends on the number of exits.

Scoring algorithm

In the previous section we assumed a situation D ∼ P(exp(Xθ + offset)). We want to find the value of θ which maximizes the penalized likelihood P (θ) = (θ) -1/2 θ P θ. We further assume that θ must satisfy a set of linear constraint Hθ = c. Then the penalized likelihood may be maximized by setting k = 0, dev 0 = ∞ and, starting from an adequate vector of initial parameters θ (0) , iterating the following steps of the so-called scoring algorithm:

k ← k + 1 d (k) = exp(Xθ (k-1) + offset) dev (k) = 2 d log d d ; STOP if |dev (k) -dev (k-1) | < tol W (k) = d (k) * I n ∆ (k) = (X W (k) X + P + H H) Ψ (k) = ∆ (k)-1 -∆ (k)-1 H (H∆ (k)-1 H ) -1 H∆ (k)-1 θ (k) = Ψ (k) X W (k) Xθ (k-1) + X (d -d (k) ) + ∆ (k)-1 H (H∆ (k)-1 H ) -1 H∆ (k)-1 c
Furthermore, by noting θ, Ŵ and Ψ the values for θ (k) , W (k) and Ψ (k) obtained at the last iteration of the scoring algorithm, the variance-covariance matrix Σ2 and the hat-matrix H are given by Σ2 = Ψ and H = X ΨX Ŵ while the effective dimension of the model df is df = tr(H) = p -q -tr( ΨP ) where p is the number of parameters (the length of θ) and q is the number of constraints (the number of rows of the constraints matrix H).

The variance of the fitted values X θ is given by var

(X θ) = diag(X ΨX ) = (X Ψ) * X 1 p
We refer the reader to Currie ( 2013) who provides a proof for the previous algorithm in the more general case of GLMs. We used the exact same notations than this author with two exceptions: our vector of parameters is called θ instead of β and the vector of estimated number of deaths at each step of the algorithm is called d (k) when he used µ (k) .

Computational aspects

In order to implement the scoring algorithm, we rely on Generalized Linear Array Models (GLAM, a framework introduced in [START_REF] Currie | Generalized linear array models with applications to multidimensional smoothing[END_REF] to speed-up computations within the scoring algorithm, especially when data has 2 or more dimensions such as our disabled mortality. The main idea behind GLAMs is to keep the matrix or array structure of the observations, offset and parameters through all the computations as opposite to standard GLMs. GLAMs allow for faster computation of linear functions Xθ (k) and X (d -d k ) as well as inner products X W (k) X and diag(X ΨX ) in the scoring algorithm or subsequent diagnosis computations.

Let us assume the design matrix may be written

X = X d ⊗ • • • ⊗ X 2 ⊗ X 1
where X i has dimensions p i × n i . We define:

• The H-transform H(X, A) of an array A by a matrix X as the array obtained by creating the matrix A * from A by flattening all dimensions except for the first one, computing A * X and reinstating the flattened dimensions. The H-transform generalize the usual matrix product in higher dimensions.

• The rotation R(A) of an array A obtained by applying a permutation to the dimensions of A so that the first dimension becomes the last.

• The rotated H-transform ρ(X, A) of an array A by a matrix X is the combination of both previous operations ρ(X, A) = R(H(X, A))

Let Θ (k) (p 1 ו••×p d ) , D (n 1 ו••×n d ) and D (k) (n 1 ו••×n d ) be the array forms of θ (k) , d and d (k) , i.e. vec Θ (k) = θ (k) , vec D = d, vec D (k) = d (k)
Then the linear functions Xθ (k) and X (d -d k )is obtained using rotated H-transforms:

•

Xθ (k) = vec ρ(X d , . . . , ρ(X 2 , ρ(X 1 , Θ (k) )) . . . ) • X (d -d k ) = vec ρ(X d , . . . , ρ(X 2 , ρ(X 1 , D -D (k) )) . . . )
Let us introduce for a matrix X (n,p) the row tensor matrix G(X) = (X ⊗ 1 p ) * (1 p ⊗ X). Then the inner product X W (k) X may be computed using

X W (k) X (p 1 ו••×p d ,p 1 ו••×p d ) ≡ ρ(G(X d ) , . . . , ρ(G(X 2 ) , ρ(G(X 1 ) , D (k) )) . . . ) (p 2 1 ו••×p 2 d )
where ≡ indicates that both expressions contain the same elements but in a different order.

In order to obtain the elements of X W (k) X in the right order:

• Computation of ρ(G(X d ) , . . . , ρ(G(X 2 ) , ρ(G(X 1 ) , D (k) )
) . . . ) should be performed

• The result must be stored in an array of dimensions p 1 × p 1 , . . . , p d × p d

• A permutation of the array dimensions must then be made with the permutation vector (1, 3, . . . , 2d-1, 2, 4, . . . , 2d) so that the odd dimensions get in front of the even dimensions

• The array must then be finally converted in a

p 1 × • • • × p d , p 1 × • • • × p d array, effectively yielding the result X W (k) X
In a similar manner, diag(X ΨX ) may be computed but in this case the reordering must be performed on the matrix Ψ in a reversed order compared to the previous case:

• Values of Ψ(p 1 ו••×p d ,p 1 ו••×p d )
must be placed in an array with dimensions (p 1 , . . . , p d , p 1 , . . . , p d )

• A permutation of the array dimensions must then be made with the permutation vector

(1, d + 1, 2, d + 2, . . . , d -1, 2d -1, d, 2d)
• Finally the array should be converted in a

p 2 1 × • • • × p 2 d array Ψr
Then we have the following relation

diag(X ΨX ) = ρ(G(X d ), . . . , ρ(G(X 2 ), ρ(G(X 1 ), Ψr )) . . . )
In the special case where the design matrix is of the form X = (X 1 :

• • • : X c ) where X i = X i,d ⊗ • • • ⊗ X i,1
and X i,j is of dimensions n i,j × p i,j , the inner product X W (k) X should be computed by blocks. For two matrices X 1 and X 2 of respective dimensions p 1 × n 1 and p 2 × n 2 , we define their row tensor matrix as

G 2 (X 1 , X 2 ) = (X 1 ⊗ 1 p 2 ) * (1 p 1 ⊗ X 2 ).
We then have:

X W (k) X =     X 1 W (k) X 1 . . . X 1 W (k) X c . . . . . . X c W (k) X 1 . . . X c W (k) X c     with X j W (k) X k (p j,1 ו••×p j,d ,p k,1 ו••×p k,d ) ≡ ρ(G 2 (X j,d , X k,d ) , . . . , ρ(G 2 (X j,2 , X k,2 ) , ρ(G 2 (X j,1 , X k,1 ) , D (k) ))) (p j,1 ×p k,1 ,...,p j,d ×p k,d )
A similar expression may easily be obtained for diag(X ΨX ).

In order to find the optimal smoothing parameters according to BIC, we use a brute-force approach and fit the model over a grid of values for the scale of the penalty. We use scales of the form 2 k for k in a reasonable range (in practice k ∈ [-5, 20] contained all our optimal parameters with a good margin). To speed-up the process, some quantities in the algorithm that does not depend on the penalty scale are computed only once, and the fitted values for the parameters for one value of the scale are used as starting values for the next scale. In practice the code required in this paper can be run in under 5 minutes on an ordinary office computer.

Autonomous mortality and incidence in LTC

We rely on data from a French private insurance portfolio containing information about more than 200,000 contributors with over 20,000 observed claims. For autonomous mortality and incidence in LTC we perform simple smoothing over the age dimension. The transition intensity is modelled by log λ x = α x where α x = B x θ x with θ x a vector of size nb x contains the scaling parameters for the B-splines and B x is the associated regression matrix. The penalty matrix is P = π x D 2,nbx D 2,nbx and there is no constraint in this case. The previously introduced scoring algorithm is used to fit the model for different values of π x . Figure 5 gives the value of π x which maximizes the BIC and has been selected as the optimal choice. Figures 6 shows the fitted transition intensities for autonomous mortality and incidence in LTC, the associated deviance residuals and extrapolation for all ages between 50 and 120. To preserve confidentiality of business-related information, the exposure in our data has been multiplied by an arbitrary coefficient, which results in a simple translation of the α x coefficients on the logarithmic scale and no impact on the results or their interpretation.

The autonomous mortality increases almost exponentially for both genders and is significantly higher for males than for females. This behaviour reminds of the general mortality with which the autonomous mortality is indeed highly correlated. As regards the incidence in LTC it also has an exponential shape for both genders but males and females incidence are much closer in this case with a higher incidence rate for males before age 85 and a slightly higher incidence rate for female after 85. The linear shape of the extrapolation is a direct consequence of our choice of an order 2 penalty. Given the linear shape of the fitted curve, we feel confident about using a linear asymptote for the extrapolation. Finally let us point out that the Mort1Dsmooth function from the MortalitySmooth [START_REF] Camarda | Mortalitysmooth: An R package for smoothing poisson counts with P-splines[END_REF] package on the statistical software R (R Core Team, 2019) could have been used directly to obtain similar results.

Mortality in LTC

While the inference of autonomous mortality and incidence in LTC is pretty straightforward, additional work is required for mortality in LTC which cannot be reasonably modelled without accounting for time spent in LTC. We nevertheless start this section by applying the same one-dimensional smoothing methodology that for other biometric laws. Figure 7 provides the selected smoothing parameters and illustrate the results of this approach by giving the fitted and extrapolated rates. Mortality in LTC is consistently higher for males than for females. However, a look at the residuals on the top left corner of Figure 16 is enough to see that there is a strong impact of duration of LTC that cannot be ignored.

Non-parametric approach

Once the duration in LTC has been identified as a relevant dimension to study the risk, a natural approach is to perform smoothing over both the age and duration in LTC dimensions. We refer the reader to [START_REF] Currie | Smoothing and forecasting mortality rates[END_REF] for more details about smoothing in 2 dimensions using P-splines. We note AxD the associated model. In this case, we have log λ x,z = A x,z with A x,z = B x,z θ x,z where θ x,z is of size nb x ×nb z and B x,z = (B z ⊗B x ). The penalty is the sum on a row penalty and a column penalty P = P x ⊗ I nbz + I nbx ⊗ P z where P x and P z are one-dimension penalties with respectives scaling constant π x and π z . The scoring algorithm is used to fit the model. The Mort2Dsmooth function from the MortalitySmooth [START_REF] Camarda | Mortalitysmooth: An R package for smoothing poisson counts with P-splines[END_REF] package could also have been used to obtain those results.

As mortality in LTC sharply decline during the first year in LTC, the assumption of constant transition intensity over one-year intervals is not realistic. We therefore use a finer level of granularity by computing exposure and number of death on a monthly basis for the first year in LTC. There is no obstacle to use the P-spine methodology with different levels of granularity in the data. However, P-splines work on the transition intensity at the middle point of the one-year intervals λ x+1/2 . The 1/2 index can be easily omitted when working with annual rates because the distance between two consecutive observations is 1 no matter what. However when mixing monthly and yearly rates, one needs to pay attention as the distance between the rates for the 11th month and for the second year is 13/24 and not 1/12.

Parametric approach

In this parametric approach, we aim at finding patterns for the role of the age and duration components by restricting the interaction between them. The idea would be to a have a simpler model in the spirit of Clayton and Schifflers (1987a) or [START_REF] Lee | Modeling and forecasting US mortality[END_REF]. We start with a very basic decomposition and then consider a second model with a non-linear interaction term.

Age Duration model

The Age-Duration (AD) model is given by log λ x,z = α x + γ z where α x = B x θ x and γ z = B z θ z , θ x and θ z being of respective sizes nb x and nb z . We have log E(D x,z ) = Xθ + log e c x,z with penalty matrix P and under constraints Hθ = c where The constraint has been introduced to make the model identifiable, otherwise an infinite number of solutions may be generated by simply translating both α x and γ z in opposite directions. Figure 12 represents the fitted coefficients and their extrapolations. The γ z plot shows a sharp decline of mortality in LTC during the first few months followed by a rather flat slope for males and a slow but steady increase for females. Regarding the extrapolation, we still use an order 2 penalty for extrapolation on the age component. As regards the duration component we prefer to use a mixed penalty for the side-effect of the horizontal asymptote. It is equivalent to making an assumption that there is some value of duration in LTC after which it no longer impacts the risk. As our data cover duration in LTC up to 8 years for males and 10 years for females, this assumption has a limited impact on the risk as only few LTC claims last more than 8 or 10 years. After 12 years the slope on γ z is flat which means there is no need to compute rates for higher durations in LTC, which proves convenient to build experience tables.

X = (1 z ⊗ B x : B z ⊗ 1 x ) P = P x 0 0 P z θ = θ x θ z H = 0 * 1 x B x 1 z B z c = 0.

Augmented Age Duration model

The duration effect γ z may be interpreted as a mixture effect linked to the presence of several groups of pathologies in the population [START_REF] Biessy | Semi-Markov modeling of the loss of autonomy among elderly people: application to long-term care insurance[END_REF] with different associated mortality levels.

In this regard, one can expect the relative weight of those groups in the population of newly disabled to vary with age at entry in LTC. To take this phenomenon into account we add an additional set of parameters β x to modulate the amplitude of the duration effect by age. We call the resulting model Augmented Age-Duration (AAD) model, defined by log λ x,z = α x + β x γ z where α x = B x θ x and γ z = B z θ z . Due to the presence of a non-linear interaction term, it can no longer be fitted by using a single GLMs. Instead, the scoring algorithm must be applied to alternatively update the following GLMs: In this case an additional constraint is required and we take the mean value of β x across the sample to be 1 in order for γ z to take similar values than in the Age-Duration model. Figure 14 represent the fitted values for the AAD model. The coefficients β x decrease with respect to age, the coefficient at age 75 being approximately 50 % higher than the coefficient at age 95 for both genders. The other parameters α x and γ z retain the same shape as in the AD model.

GLM 1 : log E(D x,z ) = [γ z ⊗ B x ] X 1 θ 1 + 1 z ⊗ α x + log e c x,z offset for GLM 1 GLM 2 : log E(D x,z ) = 1 z ⊗ B x : B z ⊗ β x X 2 θ 2 + log e c x,z offset for GLM 2 with constraints H 1 θ 1 = c 1 , H 2 θ 2 = c 2 where θ 1 = θ 1 x P 1 = P 1 x H 1 = (1 x B x ) c 1 = n x θ 2 = θ 2 x θ z P 2 = P 2 x 0 0 P z H 2 = 0 * 1 z B z 1 z B z c 2 = 0.
Regarding the extrapolation we use the same order of penalties as in the age duration model for both the α x and γ z components. For β x , we use a mixed penalty as the amplitude of the duration effect is bounded as the weights of the different pathology groups are themselves bounded between 0 and 1.

Models comparison

A comparison of the different models used in this section is provided by Figure 15 Using the AAD model over the AxD model also provides additional advantages. Firstly, empirical rates may be obtained in the AAD model while in the AxD model many cells contain less than 10 observed death, if any. Then, the AAD model allows for decomposition of mortality in LTC on several 1-dimension scales. It proves especially useful when comparing experience data to a reference mortality law. Indeed when working with the AxD model, the only natural way of positioning the experience according to the reference would be to define either a single coefficient for every cell or a full of matrix of coefficients. With the AAD model this comparison may be done on each of α x , β x and γ z , allowing for intermediate solutions. One may for example wish to use the same β x and γ z as in the reference and determine the α x parameters on the experience data.

Discussion

In this paper, we introduce a framework for modelling the Long-Term Care risk where 3 transition intensities need to be inferred: the incidence in LTC as well as the active an disabled mortality. We assume observed deaths or entries in LTC to follow a Poisson distribution with log link and specify a model for transition intensities directly, considering central exposure to risk as an offset. This gives us a convenient framework to infer each of the transition intensity, that we combine with smoothing methodology based on penalized basis splines. The disabled mortality depends not only on age but also on time already spent in the LTC state, requiring the use of 2-dimension smoothing and we provide two different approaches to tackle this issue: a direct non-parametric 2-dimension smoothing and an approach with restricted interaction terms between the age and duration in LTC component. Inside this second approach, we start by removing any interaction between the two components and then further refine the model into an Augmented Age Duration model (AAD) by allowing for a multiplicative interaction term between age and duration in LTC. The AAD model proves accurate and parsimonious and offers a clearer view of the underlying effects that take part in mortality in LTC. It thus offers more flexibility when comparing experience data to a reference table.

Compared to the approach developed in [START_REF] Tomas | Multidimensional smoothing by adaptive local kernel-weighted log-likelihood: Application to long-term care insurance[END_REF], our approach offers additional flexibility as the granularity of the data may consist of both monthly and yearly rates. The AAD model allows for easy extrapolation of the mortality rates on both the age and duration dimension, a real strength compared to non-parametric methods such asGuibert and Planchet (2018). While Biessy (2017) also provides similar advantages such as parsimony and interpretation of parameters, as a fully parametric model it comes with restrictive assumptions that are absent from the AAD model.

This approach was illustrated on the case of LTC but it may naturally be extended to other biometric risks such as Disability or Critical Illness. It could also be used for the inference of lapse rates in any life insurance portfolio. The methodology we introduce in this paper is inspired from an innovative formulation of the [START_REF] Lee | Modeling and forecasting US mortality[END_REF] model with smoothing of coefficients using P-splines [START_REF] Eilers | Flexible smoothing with B-splines and penalties[END_REF] that was developed in [START_REF] Currie | Smoothing constrained generalized linear models with an application to the lee-carter model[END_REF]. The Lee-Carter model is a stochastic model which aims at forecasting mortality. In this paper we traded the calendar year dimension for duration in LTC. It is however possible to develop a stochastic version of the introduced models by bringing back the calendar year dimension, making it a 3-dimension inference problem with smoothing on the age and time spent in LTC dimensions and no smoothing on the calendar year dimension. This stochastic model extension will be explored in further research.

Figure 2

 2 Figure 2: B-splines of degrees 0 to 3 and their derivatives

Figure 3 :Figure 4 :

 34 Figure 3: Basis of B-splines of degree 3 which consists in 5 internal knots, 6 external knots and 7 splines for observations in the range from x = 3 to x = 7.

Figure 5 :Figure 6 :

 56 Figure 5: Optimal smoothing parameters for autonomous mortality (top) and incidence in LTC (bottom) for males (left) and females (right)

Figure 7 :

 7 Figure 7: Age model; top: smoothed age coefficients for males (blue) and females (red), fitted (left) and with extrapolations (right); bottom: selected smoothing parameters

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure8: Age x Duration model for males, with minimal smoothing (left) and optimal smoothing according to BIC (right). Top: males; Bottom: females Figure8show the results of this two-dimensional smoothing approach for the selected smoothing parameters represented on Figure9. For clarity sake, the left panel of those figures represent a slightly smoothed version of the empirical rates who are very volatile as some cells contains very few observations. Mortality in LTC is very high during the first year in LTC, for males and females, especially for lower ages at entry in LTC. Residuals on the top right corner of Figure16show no obvious pattern. Finally, Figure10provides an extrapolation of the mortality surfaces for ages 50 to 120 and duration up to 15 years.

Figure 12 :

 12 Figure 11: Age Duration model: optimal smoothing parameters according to BIC

Figure 13 :

 13 Figure 13: Augmented Age Duration model: optimal smoothing parameters according to BIC

Figure 14 :

 14 Figure 14: Augmented Age Duration model: age (top) duration (bottom) and age-modulating (middle) coefficients for males (blue) and females (red). Left: fitted values; right: extrapolations

Figure 15 :

 15 Figure 15: Comparison of mortality in LTC for the 4 different models; Top: males; Bottom: females

Figure 16 :

 16 Figure 16: Residuals of mortality in LTC for the 4 different models; Top: males; Bottom: females

  whereas 16 contains the associated residuals. The AxD model obtained by 2-dimension smoothing presents the most features while the age model is clearly missing paramount duration related effects. The AD and AAD models lie in-between. With the exception of the age model, all models exhibit very similar residuals. Table1provide the values of deviance, degrees of freedom and BIC for all 4 models. Looking at deviance or BIC the age model does a very poor job at estimating mortality in LTC. The AD model is less accurate than the AxD model but compensates with a lower degree of freedom. According to BIC, it even proves a better choice than the AxD model for females. The AAD model combines a low deviance, roughly the same as the AxD model with a low degree of freedom, very close to the AD model. According to BIC or other metrics based on deviance, it therefore proves the best choice of the four models.

Table 1 :

 1 Comparison of models for mortality in LTC

	Model Gender deviance	df	BIC
	A	male	1632.58	3.05 1651.64
	AxD	male	562.63 16.15	663.65
	AD	male	619.76 10.56	685.82
	AAD	male	560.36 11.73	633.74
	A	female	1538.31	3.31 1559.32
	AxD	female	598.98 21.19	733.46
	AD	female	632.61 13.86	720.56
	AAD	female	598.33 14.74	691.88

Contact: gbiessy@scor.com.

SCOR Global Life SE -5 avenue Kléber, 75795 Paris Cedex 16, France.

2


Ben Ramzi

Mhenni

Sebastien Bourguignon

Jordan Ninin

Global Optimization for Sparse Solution of Least Squares Problems

Keywords: Sparse approximation, Cardinality constraint, Branch-and-bound, Continuous relaxation, Homotopy continuation

Finding solutions to least-squares problems with low cardinality has found many applications, including cardinality-constrained portfolio optimization, subset selection in Statistics, and many sparsity-enhancing inverse problems in signal processing. In general, this problem is NP-hard, and most works from a global optimization perspective consider a mixed integer programming (MIP) reformulation with binary variables, whose resolution is performed via branch-and-bound methods. We propose dedicated branch-and-bound algorithms for three possible formulations: cardinality-constrained and cardinality-penalized least-squares, and cardinality minimization under quadratic constraints. We show that the continuous relaxation problems involved at each node of the search tree are 1-norm-based optimization problems. A dedicated algorithm is built, based on the homotopy continuation principle, which efficiently computes the relaxed solutions for the three kinds of problems. The performance of the resulting global optimization procedure is then shown to compete with or improve over the CPLEX MIP solvers, especially for problems involving quadratic constraints. The proposed strategies are able to exactly solve some problems involving 500 to 1 000 unknowns in less than 1 000 seconds, for which CPLEX mostly fails.

Introduction

We are interested in solving optimization problems mixing a quadratic data adjustment term and a sparsity measure. Such problems arise in many application fields, among which portfolio optimization (see, e.g., [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] Cui | Convex relaxations and miqcqp reformulations for a class of cardinality-constrained portfolio selection problems[END_REF][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]), sparse regularization for inverse problems [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Elad | Sparse and Redundant Representations[END_REF][START_REF] Mendel | Optimal Seismic Deconvolution[END_REF][START_REF] Soussen | From bernoulli gaussian deconvolution to sparse signal restoration[END_REF][START_REF] Zala | High-resolution inversion of ultrasonic traces[END_REF] and compressed sensing [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF][START_REF] Elad | Compressed Sensing: Theory and Applications[END_REF], and variable or subset selection in statistics [START_REF] Efron | Least angle regression[END_REF][START_REF] Miller | Subset selection in regression[END_REF][START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. In operations research, many works addressed the cardinality-constrained problem:

P 2/0 : min x∈R n 1 2 y -Ax 2 2 subject to x 0 ≤ K,
where y ∈ R m and A ∈ R m×n with usually n > m, ||x|| 0 = Card{i|x i = 0} (which will be called the 0 -"norm" in this paper) and K ∈ N is a given cardinality that is fixed a priori. In some applications, however, one may prefer solving the error-constrained problem:

P 0/2 : min x x 0 s.t. 1 2 y -Ax 2 2 ≤ .
Such a formulation may be more relevant in signal processing and statistical problems, where parameter ≥ 0 controls the approximation level (or the prediction accuracy), that can be tuned according to prior knowledge on the data. One is then interested in finding the sparsest approximation compatible with some given noise level (or some prediction quality requirement). Finally, the penalized problem:

P 2+0 : min x 1 2 y -Ax 2 2 + µ x 0 ,
where µ > 0 trades off between approximation error and sparsity, is also encountered in the field of inverse problems, e.g., for Geophysics [START_REF] Mendel | Optimal Seismic Deconvolution[END_REF] or ultrasonic non-destructive testing [START_REF] O'brien | Recovery of a sparse spike time series by L1 norm deconvolution[END_REF][START_REF] Zala | High-resolution inversion of ultrasonic traces[END_REF]. In the Bayesian statistical framework, the 0 -norm penalization term corresponds to a Bernoulli-Gaussian prior assumption on the unknown components in x, and parameter µ then depends on the noise variance and the expected rate of non-zero values in x [START_REF] Soussen | From bernoulli gaussian deconvolution to sparse signal restoration[END_REF].

Due to the discrete nature of the 0 norm, these three problems are essentially combinatorial. Therefore, many works in signal processing and statistics proposed local optimization methods. On the one hand, we find the substitution of the 0 norm by the 1 norm, for which the three problems resort to convex, non-smooth optimization. In statistics, the 1 -norm version of problem P 2/0 is known as the Least Absolute Shrinkage and Selection Operator (LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]). On the other hand, greedy algorithms (e. g., [START_REF] Chen | Orthogonal least squares methods and their application to nonlinear system identification[END_REF][START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF][START_REF] Soussen | From bernoulli gaussian deconvolution to sparse signal restoration[END_REF]), which can similarly solve the three forms above, have been developed. In very special cases as those addressed within the compressed sensing theory [START_REF] Elad | Compressed Sensing: Theory and Applications[END_REF], it has been shown that the 1 -norm approach or greedy methods may give the solution to the 0 problem. However, in general, these approaches do not provide the global solution to the initial problem with the 0 norm [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF], and the three considered optimization problems are usually NPhard [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF].

Dedicated global optimization algorithms have been proposed for solving P 2/0 , which is a cardinality-constrained quadratic program (CCQP), in the context of sparse portfolio selection and subset selection problems. To our knowledge, Bienstock was the first to propose a specific branchand-cut algorithm for such problems (including positivity constraints) [START_REF] Bienstock | Computational study of a family of mixed-integer quadratic programming problems[END_REF]. Continuous relaxation problems involved at each iteration were solved via a specific convex quadratic programming algorithm. In [START_REF] Bertsimas | Algorithm for cardinality-constrained quadratic optimization[END_REF], Bertsimas and Shioda extended this work using Lemke's pivoting method to solve the continuous relaxation. When the matrix involved in the quadratic term is the sum of a diagonal positive matrix and a positive definite one, branch-and-bound techniques using perspective reformulation [START_REF] Frangioni | Perspective cuts for a class of convex 0-1 mixed integer programs[END_REF], Lagrangian relaxation [START_REF] Cui | Convex relaxations and miqcqp reformulations for a class of cardinality-constrained portfolio selection problems[END_REF][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] or geometric approaches [START_REF] Gao | Optimal cardinality constrained portfolio selection[END_REF], were shown to give tighter lower bounds than the continuous relaxation. In many problems of the form P 2/0 , however, variables involved in the columns of matrix A (explaining variables in statistics, dictionary atoms in sparse approximation) are often highly correlated, so that matrix A T A is ill-conditioned and such a decomposition is not possible. Much fewer works could be found that address problems P 0/2 and P 2+0 from a global optimization perspective. We remark, however, that P 2/0 and P 2+0 both involve a quadratic cost function, therefore dedicated optimization methods developed for the former may also be adapted to the latter. This is not the case for P 0/2 , for which quadratic constraints may require adapted methodology.

In [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF], the three problems were reformulated as mixed-integer programs (MIPs). Resolution was performed with the CPLEX solver. It was experimentally shown that exact resolution of difficult sparse approximation problems was possible in practice for small-size problems, but that computing times for solving P 0/2 were far higher than for the two other ones. In this paper, we propose a branch-and-bound resolution strategy for such problems. Following the works in [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], our motivation lies on the fact that sparsity-enhancing least-squares problems are very specific MIPs, that could be advantageously solved by dedicated implementation. In particular, we show that all continuous relaxation problems involved in the resolution are particular forms of 1 -norm-based problems, for which we build a dedicated continuous optimization method, based on the homotopy principle [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF][START_REF] Efron | Least angle regression[END_REF][START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF]. Our algorithm is able to solve the relaxation problems involved in any of the three problems P 2/0 , P 0/2 and P 2+0 . Therefore, these three problems can be solved exactly by our procedure, with similar computational burden.

The paper is organized as follows. In Section 2, we study the structure of continuous relaxation problems involved at each node of a branch-and-bound algorithm solving P 2/0 , P 0/2 and P 2+0 , which are reformulated as 1 -norm-based problems. In Section 3, we build a dedicated algorithm for such problems, based on the homotopy principle. In Section 4, the performance of our method is evaluated through numerical experiments, on both sparse deconvolution and subset selection problems. The discussion in Section 5 closes the paper.

Continuous relaxations within a branch-and-bound algorithm

In this section, we consider the MIP reformulations of problems P 2/0 , P 0/2 and P 2+0 , and we study/reformulate the continuous relaxation problems involved at each node of the branch-andbound strategy.

We introduce binary decision variables b i such that b i = 0 ⇔ x i = 0. We use the classical bigM formulation: assuming that solutions of interest satisfy ∀i, |x i | ≤ M for some known value M , the former logical constraint reads

-M b i ≤ x i ≤ M b i . Trivial extensions of such assumption write -M inf i ≤ x i ≤ M sup i , with M inf i , M sup i ≥ 0
, but in the following we keep |x i | ≤ M to simplify notations. The three problems can then be reformulated as the standard MIPs given in Table 1 (see for example [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF]).

Problem

MIP Reformulation

min x∈R n 1 2 y -Ax 2 2 s.t. x 0 ≤ K x ∞ ≤ M min b∈{0,1} n ,x∈R n 1 2 y -Ax 2 2 s.t. n i=1 b i ≤ K |x| ≤ M b min x∈R n x 0 s.t. 1 2 y -Ax 2 2 ≤ x ∞ ≤ M min b∈{0,1} n ,x∈R n n i=1 b i s.t. 1 2 y -Ax 2 2 ≤ |x| ≤ M b min x∈R n 1 2 y -Ax 2 2 + µ x 0 s.t. x ∞ ≤ M min b∈{0,1} n ,x∈R n 1 2 y -Ax 2 2 + µ n i=1 b i s.t.
|x| ≤ M b Table 1. Initial problems (left) and their MIP reformulations (right).

We consider a resolution strategy based on a branch-and-bound procedure, as adopted by most MIP solvers. The initial problem defines a root node. At each iteration of the algorithm, one node is selected from the list of subproblems that have not been processed yet, and a lower bound for the node is computed via the continuous relaxation of the binary variables. If this bound is greater than the current upper bound, defined by the best feasible solution found, then the subproblem is discarded. Otherwise, two children of this node are built through the addition of constraints fixing one of the relaxed variables b i to 0 and 1. The two new nodes are then added to the list. In the following, we focus on the node evaluation step by continuous relaxation.

Continuous relaxation of the root node

At the root node, no decision has been made concerning any binary variable. The continuous relaxation of binary variables in P 2/0 then reads:

P R 2/0 : min b∈[0,1] n ,x∈R n 1 2 y -Ax 2 2 s.t. n i=1 b i ≤ K |x| ≤ M b .
Proposition 1. Let P 2/1 be the following problem:

P 2/1 : min x∈R n 1 2 y -Ax 2 2 s.t. x 1 ≤ KM x ∞ ≤ M .
Then, P R 2/0 and P 2/1 have the same minimum value.

Proof. Let (b * R , x * R ) be a minimizer of P R 2/0 and let x * 1 be a minimizer of

P 2/1 . Let b = 1 M x * 1 . Then, ( b, x * 1 ) is a feasible solution of P R 2/0 , therefore y -Ax * R 2 2 ≤ y -Ax * 1 2 2 . Conversely, consider b R 2/0 = 1 M |x * R |. Then, x * R is a feasible solution of P 2/1 because x * R 1 = M b 1 ≤ KM and x * R ∞ ≤ M . Consequently, y -Ax * 1 2 2 ≤ y -Ax * R 2 2 .
We note that the same result was given in [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] for problems with positivity constraints, and in [START_REF] Bertsimas | Best subset selection via a modern optimization lens[END_REF] in our case. A similar result holds for P R 0/2 , the continuous relaxation of P 0/2 :

P R 0/2 : min b∈[0,1] n ,x∈R n n i=1 b i s.t. 1 2 y -Ax 2 2 ≤ |x| ≤ M b
Proposition 2. Let P 1/2 be the following problem:

P 1/2 : min x∈R n x 1 s.t. 1 2 y -Ax 2 2 ≤ x ∞ ≤ M .
Then, P R 0/2 and P 1/2 have the same minimum value.

Proof. Let (b * R , x * R ) an optimal solution of P R 0/2 . Suppose that |x * R i | < M b * R i for some compo- nent i. Consider b = 1 M |x * R |, such that b i < b * R i . Then, (b , x * R ) is a feasible solution for P R 0/2 with n i=1 b i < n i=1 b * R i , which contradicts the definition of (b * R , x * R ). Indeed, |x * R | = M b * R and x * R 1 = M n i=1 b * R i . Thus (b * R , x * R ) is an optimal solution of P 1/2 .
Finally, consider the continuous relaxation of binary variables in the penalized problem P 2+0 as follows:

P R 2+0 : min b∈[0,1] n ,x∈R n 1 2 y -Ax 2 2 + µ n i=1 b i s.t.|x| ≤ M b .
Proposition 3. Let P 2+1 be the following problem:

P 2+1 : min x∈R n 1 2 y -Ax 2 2 + µ M x 1 s.t. x ∞ ≤ M.
Then, P R 2+0 and P 2+1 have the same minimum value.

Proof. The proof is similar to that of Proposition 2.

The three continuously relaxed problems can therefore be reformulated without binary variables, and correspond to 1 optimization problems with box constraints.

Continuous relaxation in the branch-and-bound algorithm

We now consider a given node in the branch-and-bound algorithm and the corresponding continuous relaxation sub-problem. Let S 0 (respectively, S 1 ) denote the index set of binary variables that are set to 0 (respectively, to 1), and let S collect all remaining indices:

∀i ∈ S 1 , b i = 1 and |x i | ≤ M ∀i ∈ S 0 , b i = 0 and x i = 0 ∀i ∈ S, b i ∈ [0, 1] and |x i | ≤ M b i
Let A S denote the sub-matrix formed by all columns of matrix A indexed by S. Similarly, u S denotes the corresponding sub-vector of u.

For problem P 2/0 , the continuous relaxation of variables bS in the corresponding sub-problem reduces to:

Q R 2/0 : min x S 1 ∈ R n1 bS ∈ [0, 1] n xS ∈ R n 1 2 y -A S 1 x S 1 -ASxS 2 2 s.t.      i∈ S b i ≤ K -n 1 |xS| ≤ M bS x S 1 ∞ ≤ M ,
where n 1 and n denote the size of S 1 and S. Then, similarly to the developments in Section 2.1, one can show that Q R 2/0 and Q 2/1 have the same minimum value:

Q 2/1 : min x S 1 ,xS 1 2 y -A S 1 x S 1 -ASxS 2 2 s.t.      xS 1 ≤ M (K -n 1 ) xS ∞ ≤ M x S 1 ∞ ≤ M .
Applying a similar reasoning to the two other formulations, we finally obtain the equivalent problems summarized in Table 2.

Two important comments arise from these results:

• Whatever the formulation (constrained or penalized), all continuously relaxed subproblems involved in the evaluation of the nodes in the branch-and-bound algorithm can be reformulated without binary variables. • They all reduce to optimization problems mixing a least-squares function, 1 -norm terms involving only a part of the variables, and box constraints.

In Section 3, we build a dedicated algorithm which solves the three problems Q 2/1 , Q 1/2 and Q 2+1 .

A dedicated homotopy continuation algorithm for relaxed problems

Optimization involving quadratic misfits and 1 -norm terms has been a very active field of research in the past ten years. Many dedicated convex, non-smooth, optimization algorithms have been

Continuous relaxation problem

Equivalent

1 -norm-based problem Q R 2/0 : min b∈[0,1] n ,x∈R n 1 2 y -Ax 2 2 s.t. n i=1 b i ≤ K |x| ≤ M b b S 1 = 1 b S 0 = 0 Q 2/1 : min x S 1 ,xS 1 2 y -A S 1 x S 1 -ASxS s.t. xS 1 ≤ M (K -n 1 ) xS ∞ ≤ M x S 1 ∞ ≤ M Q R 0/2 : min b∈[0,1] n ,x∈R n n i=1 b i s.t. 1 2 y -Ax 2 2 ≤ |x| ≤ M b b S 1 = 1 b S 0 = 0 Q 1/2 : min x S 1 ,xS 1 M xS 1 + n 1 s.t. 1 2 y -A S 1 x S 1 -ASxS ≤ xS ∞ ≤ M x S 1 ∞ ≤ M Q R 2+0 : min b∈[0,1] n ,x∈R n 1 2 y -Ax 2 2 + µ n i=1 b i s.t. |x| ≤ M b b S 1 = 1 b S 0 = 0 Q 2+1 : min x S 1 ,xS 1 2 y -A S 1 x S 1 -ASxS 2 2 + µ M xS 1 + µ n 1 s.t. xS ∞ ≤ M x S 1 ∞ ≤ M Table 2.
Continuous relaxation problems at any node in the branch-and-bound procedure (left), and equivalent problem without binary variables involving the 1 norm (right), for the three considered formulations. developed (see for example [START_REF] Elad | Compressed Sensing: Theory and Applications[END_REF][START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF] and references therein), for solving problems:

P τ 1 : min x 1 2 y -Ax 2 2 s. t. x 1 ≤ τ ; P 1 : min x x 1 s. t. 1 2 y -Ax 2 2 ≤ ; P µ 1 : min x 1 2 y -Ax 2 2 + µ x 1 .
Most works addressed the penalized form, which resorts to unconstrained optimization. Since the two objectives are convex, the three problems are equivalent, that is, for any τ ≥ 0, there exists µ (τ ) ≥ 0 such that P τ 1 and P µ (τ )

1 have the same solution, and reciprocally. Similarly, for any ≥ 0, there exists µ ( ) ≥ 0 such that P 1 and P µ (τ ) 1 are equivalent. However, in general, there is no explicit mapping between the three parameters.

In this section, we propose to solve the three problems Q 2/1 , Q 1/2 and Q 2+1 by a generalization of the homotopy continuation principle [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF][START_REF] Efron | Least angle regression[END_REF][START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF]. This choice is motivated by the following reasons:

• First and foremost, the same algorithm can solve the three problems with the same computational efficiency. • Then, it can naturally incorporate specificities of the problems

Q 2/1 , Q 1/2 and Q 2+1 , such
as box constraints and partial 1 -norm terms. • In addition, it is an exact algorithm, for which the solution is obtained in a finite number of iterations. • Last, its efficiency increases with the sparsity level in the sought solution.

The homotopy continuation method considers the penalized form P µ 1 , and exploits the fact that the solution path is piecewise linear as a function of µ. Starting from µ (0) = A T y ∞ (such that the solution is identically zero ∀µ > µ (0) ), it iteratively computes all solutions by continuously decreasing parameter µ until the target value is reached. Therefore, it can similarly solve problems P τ 1 or P 1 , by stopping when the corresponding value of τ or is reached, respectively. Figure 1 shows a typical solution path (left part) and the corresponding evolution on the Pareto front (right part).

µ x * 1 x * 2 x * 3 x * 4 x * µ (0) µ (1) µ (2) µ (3) µ (4) µ * 1 2 y -Ax * 2 x * 1 τ µ (0) µ (1)
µ (2) µ (3) µ (4) µ * Solution path for P µ 1 In the following, we generalize the homotopy method to the class of problems Q 2/1 , Q 1/2 and Q 2+1 , where the 1 norm only operates on some of the variables and box constraints are included. Note that the homotopy method with box constraints was recently proposed in [START_REF] Liang | Homotopy algorithm for box-constrained LASSO and its convergence[END_REF], which also established convergence proofs. Including partial 1 -norm penalization mainly impacts initialization and additional tests to be performed at each iteration. We first derive optimality conditions in § 3.1. Then, the algorithm is built in § 3.2. Implementation details are discussed in § 3.4.

Corresponding Pareto curve

Optimality conditions

We first focus on problem Q 2+1 in Table 2, that is, the relaxed problem involved in the cardinalitypenalized form. We remove from the problem the variables x S 0 which are fixed to zeros and consider equivalently the optimization problem of the form:

min x F (x) := J(x) + λh(x) s. t. g i (x) ≤ 0 ∀i = 1, . . . , n, with J(x) := 1 2 y -ASxS -A S 1 x S 1 2 2 h(x) := xS 1 g i (x) := |x i | -M . (1)
The Lagrangian function is:

L(x, π) = J(x) + λh(x) + n i=1 π i g i (x),
with π ∈ R n the vector of Lagrange multipliers associated with the bound constraints g i (x) ≤ 0. The function J is differentiable with ∇J(x) = -A T (y -Ax). The subdifferentials of h(x) = xS 1 and of g i (x) = |x i | -M are respectively:

∂h(x) =    z ∈ R n z i = 0 if i ∈ S 1 ∪ S 0 z i = sgn(x i ) if i ∈ S and x i = 0 z i ∈ [-1, 1] if i ∈ S and x i = 0    (2)
and

∂g i (x) =    z ∈ R n z j = 0 for j = i z i = sgn(x i ) if x i = 0 z i ∈ [-1, 1] if x i = 0    . (3)
The vector x is optimal if and only if there exists π such that (x , π) satisfies the Karush-Kuhn-Tucker optimality conditions applied to continuous convex non-differentiable functions (see for example [START_REF] Rockafellar | Convex Analysis[END_REF]):

(P KKT )                0 ∈ -A T (y -Ax ) + λ∂h(x ) + ∂ n i=1 π i g i (x ) (4)
g i (x ) ≤ 0 ∀i = 1, . . . , n (5)
π i ≥ 0 ∀i = 1, . . . , n (6)
π i g i (x ) = 0 ∀i = 1, . . . , n. (7)
Particular points are those which activate the bound constraint (x i = ±M for i ∈ S ∪ S 1 ) or nondifferentiability points (x i = 0 for i ∈ S). We now split the variable indices into the five possible cases.

-Case 1. Let S0 = {i ∈ S | |x i | = 0}. From equation [START_REF] Chen | Orthogonal least squares methods and their application to nonlinear system identification[END_REF], πS 0 = 0 and using equations ( 2) and ( 4), optimality conditions of x S0 become:

|A T S0 (y -A S 1 x S 1 -ASx S)| < λ. ( 8
) -Case 2. Let Sin = {i ∈ S | 0 < |x i | < M }.
From equation [START_REF] Chen | Orthogonal least squares methods and their application to nonlinear system identification[END_REF], πS in = 0 and using equations ( 2) and ( 4), optimality conditions of x Sin become:

-A T Sin (y -A S 1 x S 1 -ASx S) + λsgn(x Sin ) = 0. (9)
-Case 3. Let S = {i ∈ S | |x i | = M }.
From equation ( 6), πS ≥ 0 and using equations ( 3) and ( 4), optimality conditions of x S become:

-A T S (y -A S 1 x S 1 -ASx S) + λsgn(x S ) + πS sgn(x S ) = 0, πS ≥ 0, ( 10
)
where denotes the Hadamard (entrywise) product.

-Case 4. Let S 1 in = {i ∈ S 1 | 0 ≤ |x i | < M }.
From equation ( 7), π S 1 in = 0 and using equations ( 2) and (4), optimality conditions of x S 1 in become:

-A T S 1 in (y -A S 1 x S 1 -ASx S) = 0. ( 11
) -Case 5. Let S 1 = {i ∈ S 1 | |x i | = M }.
From equation ( 6), π S 1 ≥ 0 and using equations ( 3) and ( 4), optimality conditions of x S 1 become:

-A T S 1 (y -A S 1 x S 1 -ASx S) + π S 1 sgn(x S 1 ) = 0, π S 1 ≥ 0. (12)
Let us remark that equations ( 9), ( 10) and ( 12) concern non-zero variables, therefore the corresponding sign function is well-defined.

Homotopy continuation algorithm

We now build the homotopy algorithm that solves the problem [START_REF] Bertsimas | Best subset selection via a modern optimization lens[END_REF]. Let r = y -AS x S -A S 1 x S 1 , where each component in xS and x S 1 equals ±M . Equations ( 9) and ( 11) are linear systems in x Sin and x S 1 in , whose the solution is:

   x Sin = A T Sin (I -P S 1 in )AS in -1 A T Sin (I -P S 1 in )r -λ sgn(x Sin ) , (13)
x

S 1 in = (A T S 1 in A S 1 in ) -1 (A T S 1 in r -A T S 1 in AS in x Sin ), (14)
where

P S 1 in = A S 1 in (A T S 1 in A S 1 in ) -1 A T S 1
in and I is the identity matrix of appropriate size. These equations show that, in a given configuration of the index sets { Sin , S 1 in , S , S 1 , S0 } (which we call the support configuration), the solution of the problem (1) is linear in λ. Indeed, variables in S and S 1 are fixed to ±M , and variables in S0 are zero. The homotopy method then constructs the solution path (the set of all solutions as a function of λ) by identifying iteratively the different breakpoints that lead to changes in the support configuration. These breakpoints will occur at specific values of λ, for which (at least) one of the conditions in equations ( 8)-( 12) is violated. The algorithm works as follows:

(1) First, it is clear that as λ → +∞, 1 -norm-penalized variables xS are zero. In that case, other variables x S 1 are found by solving the least-squares problem: min

-M ≤x S 1 ≤M 1 2 y -A S 1 x S 1 2 2 .
We note x (0) the vector defined by

x (0) S 1 := argmin -M ≤x S 1 ≤M 1 2 y -A S 1 x S 1 2 2 , (15) x (0)
S := 0.

Equation [START_REF] Cui | Convex relaxations and miqcqp reformulations for a class of cardinality-constrained portfolio selection problems[END_REF] shows that x (0) is the solution of the problem (1) as long as λ ≥ λ (0) , with:

λ (0) := A T S (y -A S 1 x (0) S 1 ) ∞ . (17)
(2) As λ decreases below λ (0) , indices j ∈ S such that |a T j (y

-A T S 1 x (0)
S 1 )| = λ (0) leave S0 to form the new subset Sin . This new support configuration remains valid for any λ ∈ [λ (1) , λ (0) ], where λ (1) defines the next breakpoint, etc. A monotonically decreasing sequence {λ (k) } k is built iteratively, by testing all possible changes that can occur to the support configuration, and selecting the one(s) corresponding to the smallest decrease in λ. Then, the support configuration is updated, and a new breakpoint in λ is searched. Since the solution path is piecewise linear as a function of λ, the solution x (k) at the k-th breakpoint reads:

x (k) = x (k-1) + γ (k) d (k) (18) and λ (k) = λ (k-1) -γ (k) , (19)
where d (k) represents the vector of slope changes and γ (k) > 0 represents the length of interval [λ (k) , λ (k-1) ]. From equations ( 13) and ( 14), the direction d (k) is obtained by:

         d (k) Sin = (A T Sin (I -P S 1 in )AS in ) -1 sgn(x (k-1) Sin ) (20)
d (k) S 1 in = -(A T S 1 in A S 1 in ) -1 A T S 1 in AS in d (k) Sin (21) d (k) i = 0 ∀i / ∈ { Sin ∪ S 1 in }, (22)
where the last equality concerns variables that are fixed to zero or to ±M . The step size γ (k) is obtained as the smallest positive value γ > 0 such that x (k-1) + γd (k) reaches a new breakpoint. Five different cases can occur, which are detailed hereafter. We introduce the following notation, where exponentiation by k is removed to simplify it.

c = A T (y -ASx (k-1) S -A S 1 x (k-1) S 1
) and u = A T (AS in d

(k) Sin + A S 1 d (k) S 1
).

(a) A component with index ∈ S0 becomes nonzero when equality in equation ( 8) is reached. Inserting equations ( 13) and ( 14) into equation ( 8), one can show that x become positive (respectively, negative) if:

γ = λ (k-1) + c 1 -u respectively, when γ = -λ (k-1) + c -1 -u . (23)
(b) A component with index ∈ Sin becomes zero. From equation [START_REF] Li | Optimal lot solution to cardinality constrained mean-variance formulation for portfolio selection[END_REF], this occur if:

γ = -x (k-1) d (k) . (24)
(c) A component with index ∈ Sin or S 1 in yields the bound M or -M , depending on its current sign. From [START_REF] Li | Optimal lot solution to cardinality constrained mean-variance formulation for portfolio selection[END_REF], this occur if:

γ = M sgn(x (k-1) ) -x (k-1) d (k) . ( 25
)
(d) The bound constraint for some component with index ∈ S becomes inactive. This may occur when the corresponding Lagrange multiplier π = 0 in equation [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF], which yields:

γ = sgn(x (k-1) )λ (k-1) -c sgn(x (k-1) ) -u . ( 26
)
(e) The bound constraint for some component with index ∈ S 1 becomes inactive. This may occur when the corresponding Lagrange multiplier π = 0 in equation ( 12), which yields:

γ = -c u . ( 27
)
The shortest step size γ (k) is then defined as the shortest positive step among all possible ones, defined by equations ( 23)- [START_REF] Shaw | Lagrangian relaxation procedure for cardinality-constrained portfolio optimization[END_REF]. In theory, γ (k) may be obtained by several conditions above simultaneously; should this happen, the support configuration is updated correspondingly.

(3) The algorithm stops when the target λ, say λ , is reached, that is, after iteration k such that λ ∈ [λ (k) , λ (k-1) ]. Then, the optimal solution is found by:

x = x (k-1) + γ d (k) , ( 28
)
with γ = λ -λ (k) .
The homotopy algorithm is summarized in Algorithm 1. Figure 2 shows the solution path for a toy example with 5 variables: S = {1, 2, 3} and S 1 = {4, 5}.

Algorithm 1 Homotopy algorithm for solving the problem Q 2+1 in Table 2, reformulated as the problem (1) with λ = λ .

1: Set k = 0. Initialize λ (0) , x (0) by equations ( 15)-( 17) 2: while not stop condition do 3:

k ← k + 1 4:
Update d (k) by equations ( 20)- [START_REF] Miller | Subset selection in regression[END_REF] 5:

Determine the step size γ (k) as the smallest positive value among all cases of equations ( 23)- [START_REF] Shaw | Lagrangian relaxation procedure for cardinality-constrained portfolio optimization[END_REF] 6:

Compute accordingly (x (k) , λ (k) ) by equations ( 18)- [START_REF] Liang | Homotopy algorithm for box-constrained LASSO and its convergence[END_REF] 7:

Update index sets { Sin , S 1 in , S , S 1 , S0 }

If λ (k) ≤ λ , then stop. 9: return x

3.3. Solutions to constrained problems Q 2/1 and Q 1/2
As λ is continuously decreased, the 1 norm of the penalized variables x S 1 continuously increases and the least squares function

1 2 y -ASx S -A S 1 x S 1 2
2 continuously decreases. Therefore, the homotopy method can also solve the constrained problems Q 2/1 and Q 1/2 in Table 2. More precisely:

• For Q 2/1 , iterations stop at the first breakpoint when the 1 norm of the penalized variables

x (k) S 1 exceeds the threshold value τ := M (K -n 1 ). Then, in the corresponding interval [λ (k) , λ (k-1)
], the solution is given by equation [START_REF] Soussen | From bernoulli gaussian deconvolution to sparse signal restoration[END_REF]. By construction, there is no sign change between x (k-1) and x , and the value of γ is obtained when x 1 = τ , which happens if:

γ = τ -x (k-1) 1 sgn(x (k-1) ) T d (k) . M -M λ x * x * 1 x * 2 x * 3 x * 4 x * 5 λ (0) λ (1) λ (2) λ (3) λ (4) λ (5) λ (6) λ (7) λ * γ (0) γ (1) γ (2) γ (3) γ (4) γ (5) γ (6) γ Figure 2.
Example of solution path x giving the solutions of the problem (1) as a function of λ, with 5 variables: S = {1, 2, 3} and S1 = {4, 5}. Circles represent the events that cause a change in the support configuration. Vertical dotted lines represent the breakpoints.

• Similarly, for Q 1/2 , iterations stop at the first breakpoint such that 1 2 y -ASx [START_REF] Soussen | From bernoulli gaussian deconvolution to sparse signal restoration[END_REF] in the least squares expression, the value of γ such that

(k) S - A S 1 x (k) S 1 2 2 ≤ . Substituting equation
1 2 y -ASx * S -A S 1 x * S 1
2 2 = can be found by solving a (scalar) quadratic equation, whose solution is:

γ = r T v -(r T v) 2 -v T v(r T r -2 ) v T v , with r = y -ASx (k-1) S -A S 1 x (k-1) S 1 and v = AS in d (k) Sin + A S 1 d (k) S 1 .

Implementation and practical issues

Some practical remarks concerning the numerical implementation of the homotopy algorithm 1 are detailed in this section. First, each iteration mostly consists of solving linear systems of equations ( 20) and ( 21), whose size respectively correspond to the current number of variables in Sin and S 1 in . Since the support configuration only slightly changes between two breakpoints 1 , the matrix inverses of equations ( 20) and ( 21) can be computed recursively by performing rank-one updates. In our simulations, using the block matrix inversion formulas appeared to be the most efficient strategy.

We also remark that, for a given non-zero component, all computations of equations ( 20)-( 21) are not necesarry. In particular, if x (k-1) > 0 and d (k) > 0 (or respectively, d (k) < 0), then the only possible breakpoint that can occur for component is when x reaches the upper bound M , i.e. the case 2c in §3.2 (or respectively, 0, i.e. the case 2b). A similar reasoning is applied to negative components if x (k-1) < 0.

We conclude this section with an important remark concerning the use of the homotopy algorithm in the branch-and-bound resolution of the error-constrained problem P 0/2 . At each node, the continuous relaxation brought by the solution to problems such as Q 1/2 in Table 2 provides a lower bound on the global optimum value of P 0/2 , which is compared to the best known feasible solution. Since the objective function in Q 1/2 increases within the homotopy iterations, the homotopy algorithm is stopped as soon as the 1 norm of the current iterate exceeds the global upper bound. In this case, the node of the branch-and-bound tree is pruned.

Experimental results

We now insert the continuous relaxation algorithms built in Section 3 into dedicated branch-andbound procedures for the resolution of the three problems P 2/0 , P 0/2 , and P 2+0 . We use depth-first search and we choose to branch on the continuous variable x k with the maximum absolute value in the solution of the relaxed problem. We branch up first, that is, we first explore the branch corresponding to the decision x k = 0. This strategy is well adapted for problem P 2/0 , in which the depth limit is imposed by the cardinality constraint. It also allows one to quickly find feasible solutions for problems P 2+0 and P 0/2 with limited depth search, since we know in advance that their solutions are sparse. For problem P 2+0 , we also exploit the fact that the initialization step in the homotopy algorithm at each node evaluation (see § 3.2) produces a feasible solution, which is used in order to improve the best solution found.

In this section, we compare our homotopy-based branch-and-bound implementation, named B&B R-HOM , to:

• the same branch-and-bound exploration strategy, where the continuous relaxation problems All methods are implemented in C++ and executed on a UNIX machine equipped with 31.1Go RAM and four Intel Core i7-6600U central processing units (CPUs) clocked at 2.60 GHz. Computations are restricted to one core in order to focus on the algorithm performance, disabling parallelization capacities. All experiments are run by using the default CPLEX settings, and the maximum CPU time allowed was set to 1000 s.

Q R 2/0 , Q R 0/2 and Q R 2+0 in
In § 4.1, we study the performance of our algorithm on simulated sparse deconvolution problems typically encountered in signal processing. Then, simulated subset selection problems with random entries are considered in § 4.2.

Sparse deconvolution problems

Sparse deconvolution is a classical signal processing problem [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF]Chapter 5], which aims to estimate a sparse sequence x from filtered and noisy observations. Collecting sampled data in vector y and the unknown sparse sequence in vector x, we obtain a model of the form y = Ax + n, where A is a convolution matrix and n is a random term representing noise and model errors. Deconvolution is an ill-posed inverse problem [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF], whose resolution with sparsity-enhancing regularization is often addressed through suboptimal 1 -norm-based or greedy methods. It was shown in [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF] that exact resolution of 0 -norm problems (with the CPLEX MIP solver) achieves better solutions, obviously with higher computing times. From a practical point of view, the three problems P 2+0 , P 2/0 and P 0/2 are of interest, depending on the available information for tuning their respective parameter. We consider the problem instances2 proposed in [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF], with y ∈ R 120 , A ∈ R 100×120 , noise samples are independent, identically and normally distributed with signal-to-noise ratio SNR = 10 dB and M = 1.1 A T y ∞ . Columns of A have been normalized. The cardinality varies from K = 5 to K = 9 for P 2/0 . For problems P 0/2 and P 2+0 , the respective parameters and λ are tuned from statistical rules accounting for the noise level and the sparsity degree (see [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF] for details): if n denotes the noise vector, is tuned such that the probability P ( n 2 ≤ ) = 95 %, and λ = 2σ 2 log(1/ρ -1), where σ 2 is the noise variance and ρ = K/n.

Computational results are reported in Table 3. 3. Computational results for simulated sparse deconvolution problems. For each algorithm, columns "Time", "Nodes", and "F" respectively give the average computing time in seconds, the average number of explored nodes, and the number of problems that could not be solved in 1000s. 50 instances of each problem were considered. Averages are performed only over instances which could be solved in less than 1000s by the three algorithms.

For all problems, B&B R-HOM and B&B R-Cplex explore the same number of nodes, since they use the same branch-and-bound strategy (slight differences are observed, though, that may be due to numerical issues). However, B&B R-HOM requires significantly less execution time: about a factor between 3.6 and 9 for P 2/0 and P 2+0 , and about a factor between 20 and 30 for P 0/2 . It also successfully solves more instances in 1000s.

Unsurprisingly, MIP Cplex performs better than B&B R-HOM on P 2/0 and P 2+0 : it runs generally 2 to 3 times faster (up to 5.3 on the most difficult P 2+0 problem with K = 9), and both the resolution strategy (number of nodes) and the continuous relaxation efficiency (computing time per node) are favorable to MIP Cplex . This was expected, since the CPLEX MIP solver benefits from many additional developments that are not implemented in our branch-and-bound algorithm, such as cutting planes, heuristics, warm start, etc.. However, our algorithm outperforms MIP Cplex on all instances of P 0/2 , both in terms of number of nodes (2 to 3 times less) and CPU time. In particular, the computing time per node is reduced efficiency (computing time per node) are favorable to B&B R-HOM .

We also remark that the CPLEX MIP solver performs very badly on most P 0/2 instances. In particular, the number of explored nodes in MIP Cplex is much bigger than with the proposed branchand-bound strategy. The computing time per node with MIP Cplex is also much higher than with B&B R-HOM , confirming the improvement brought by the homotopy-based relaxation computations.

Performance profiles obtained over the 900 considered instances (50 instances for each of the 18 problems defined in Table 3) are shown in right part of Figure 3, from which we have the following conclusions:

• The three algorithms generally solve more instances in 1000 s for (m, n) = (500, 1000) than for (m, n) = (250, 500). This somehow counterintuitive result can still be explained by the compressed sensing theory [START_REF] Elad | Compressed Sensing: Theory and Applications[END_REF]: as the dimension increases, the columns of A, composed of random numbers, become less correlated, making the problems easier. • The B&B R-HOM algorithm clearly outperforms B&B R-Cplex and MIP Cplex . It is the fastest and most efficient in all solved problems.

Conclusion

We have investigated exact optimization algorithms for least squares problems with sparsity constraint. By studying the problem properties, we proposed a dedicated branch-and-bound procedure, which binary variables that are usually used for MIP reformulation. In particular, an algorithm was built for solving the continuous relaxation problems involved at any node of the search tree, which were recast as continuous non-smooth optimization problems, involving the 1 norm. Inspired by the homotopy principle, this algorithm can be applied with similar efficiency for the three addressed formulations (cardinality-constrained and cardinality-penalized least-squares, and cardinality minimization under quadratic constraints). In particular, the resulting strategy was shown to efficiently solve the quadratically constrained formulation, which is of major interest in many applications, whereas the CPLEX MIP solver resulted to be much less efficient. The method is even more efficient as the solution searched is sparse. The computational burden of the homotopy algorithm increases almost linearly with the cardinality of the relaxed solutions. In cases where the cardinality increases, other 1 -norm-based methods may be investigated, which could benefit from warm starting when evaluating the continuous relaxation at two close nodes. Lagrangian relaxation was also shown to achieve tighter relaxation than the continuous one, for such problems where the Gram matrix A T A admits a specific decomposition. Studying dedicated algorithms in this case, in particular for the two less familiar formulations P 0/2 and P 2+0 , would also deserve further works. Finally, the three formulations could be tackled jointly from a multiobjective optimization perspective, involving dedicated branch-and-bound techniques [START_REF] Ehrgott | Multicriteria Optimization[END_REF].

Figure 1 . 2 2 , x 1 )

 121 Figure 1. Homotopy method: example of solution path x (µ) = arg minx 1 2 y -Ax 2 2 + µ x 1 as a function of µ (left), and corresponding set ( 1 2 y -Ax 2 2 , x 1 ) as a function of µ (right).

  Table 2 are solved with the CPLEX 12.8 quadratic programming solver (B&B R-Cplex ); • the resolution of the MIP problems in Table 1 with the CPLEX 12.8 MIP solver (MIP Cplex ).

Table

	Problem		Branch-and-bound			MIP solver CPLEX
			B&B R-HOM		B&B R-Cplex		CPLEX 12.8
			Time	Nodes	F	Time	Nodes	F	Time	Nodes	F
			(s)	(10 3 )		(s)	(10 3 )		(s)	(10 3 )
		K = 5	1.9	1.28	0	7.7	1.28	0	3.0	1.71	0
	P 2/0	K = 7	26.1	17.89	0	141.9	17.89	2	16.6	21.51	0
		K = 9 124.3	57.37	16 448.1	57.46	30	53.8	72.04	6
		K = 5	3.9	2.01	0	37.7	2.02	0	3.2	1.98	0
	P 2+0 K = 7	26.3	10.20	0	187.9	10.22	7	7.4	9.61	0
		K = 9	91.1	31.80	14 470.7	31.87	28	17.3	23.74	2
		K = 5	3.4	1.98	0	102.4	2.05	0	25.7	6.71	0
	P 0/2	K = 7	9.4	8.26	0	274.4	8.25	13	52.7	20.48	2
		K = 9	21.3	13.14	7	485.2	13.15	35 111.3	28.36	17

Although several support configuration changes may happen simultaneously, in practice, Sin and S 1 in are usually modified by at most one element at each breakpoint, corresponding to the activation of one condition among equations (23)-[START_REF] Shaw | Lagrangian relaxation procedure for cardinality-constrained portfolio optimization[END_REF].

Data are available online at pagesperso.ls2n.fr/∼bourguignon-s/download MIP.html

by more than 2 with the homotopy algorithm. As a result, B&B R-HOM runs 4 to 8 times faster than MIP Cplex on these instances, even with a rather basic branch-and-bound strategy.

We also note that the relative efficiency of B&B R-HOM compared to other tested algorithms increases as K decreases. Indeed, when the global optimum is sparse, the solutions to the continuous relaxation problems are sparse, too. Therefore, the homotopy algorithm becomes more efficient. Indeed, the computing time increases with the number of non-zero components in the solution of the relaxed problem.

In Figure 3, the performance profiles, defined by Dolan and Moré in [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF], are visual tools to benchmark algorithms. The fraction of solved problems is represented as a function of the performance ratio: for each instance, the ratio between the computing time of each algorithm and the minimum computing time is computed. Then, the performance profiles draw the cumulative distribution functions for such ratios.

Fraction of 450 problems Performance profiles obtained over the 450 considered instances (50 instances for each of the 9 problems defined in Table 3) are shown in the left part of Figure 3, from which we have the following conclusions:

• The performance of B&B R-HOM is far above that of B&B R-Cplex , confirming the results in Table 3. • MIP Cplex solves more problems in 1000 s and is globally the fastest when considering all instances (the green curve achieves the higher value on the top-right corner). • Nevertheless, the performance of B&B R-HOM is above that of MIP Cplex in most part of the profile. For example, B&B R-HOM requires 6.7 times more time than the best algorithm for solving 90% of the problems, whereas MIP Cplex requires 12.5 times more time (see the horizontal lines in Figure 3). Similarly, in 5 times the best algorithm time, B&B R-HOM solves 86% of problems, whereas MIP Cplex only solves 76% of them (see the vertical line in Figure 3).

Simulated Subset selection problems

We now study algorithmic performance on artificial subset selection problems with random entries (see [START_REF] Bertsimas | Algorithm for cardinality-constrained quadratic optimization[END_REF][START_REF] Zheng | Successive convex approximations to cardinality-constrained convex programs: a piecewise-linear dc approach[END_REF] for similar simulations). We consider a dictionary A ∈ R m×n , with n = 2 × m, composed of independent and identically distributed random entries, which are drawn from a centered, unit-variance, normal distribution (columns have been normalized). Sparse vectors are generated with random non-zero positions, and corresponding amplitudes are drawn as u + sgn(u), where u follows a centered, unit variance, Gaussian distribution. Noise is added such that SNR = 10 dB. The problem cardinality varies from K = 5 to K = 15. Parameters M and formulation-specific parameters K, and µ are tuned similarly as in § 4.1. Note that the random nature of the dictionary makes problems easier to solve, which allows us to increase the problem dimension: we consider here n = 500 and n = 1000. Indeed, the compressed sensing theory (see for example [START_REF] Elad | Compressed Sensing: Theory and Applications[END_REF]) shows that, if the correlation between columns in A is sufficiently low, then, with high probability, the 1 -norm relaxation may solve the original 0 -norm problem.

Computational results are reported in Table 4. As in § 4.1, B&B R-HOM and B&B R-Cplex explore the same number of nodes. However, B&B R-HOM requires significantly less execution time: about a factor of 10 to 100 for P 2/0 and P 2+0 , and about a factor of 100 to 230 for P 0/2 . It also successfully solves more instances in less than 1 000s. Compared to MIP Cplex , B&B R-HOM performs better on most problems in terms of execution time: it runs generally more than 10 times faster (up to 276 times on the most difficult problems P 0/2 with K = 5 and n = 1000). For all problems, both the number of nodes and the continuous relaxation