extfieldsofstudy:
Mathematics
Computer Science

s2fieldsofstudy:
Computer Science
Mathematics

text:
Dimension free non-asymptotic bounds on the accuracy of high dimensional Laplace approximation

This note attempts to revisit the classical results on Laplace approximation in a modern non-asymptotic and dimension free form. Such an extension is motivated by applications to high dimensional statistical and optimization problems. The established results provide explicit non-asymptotic bounds on the quality of a Gaussian approximation of the posterior distribution in total variation distance in terms of the so called \emph{effective dimension} \( p_G \). This value is defined as interplay between information contained in the data and in the prior distribution. In the contrary to prominent Bernstein - von Mises results, the impact of the prior is not negligible and it allows to keep the effective dimension small or moderate even if the true parameter dimension is huge or infinite. We also address the issue of using a Gaussian approximation with inexact parameters with the focus on replacing the Maximum a Posteriori (MAP) value by the posterior mean and design the algorithm of Bayesian optimization based on Laplace iterations. The results are specified to the case of nonlinear inverse problem.

Introduction
High dimensional Laplace approximation has recently gained an increasing attention in connection with Bayesian inference for complicated nonlinear parametric models such as nonlinear inverse problems and Deep Neuronal Networks. Laplace approximation is obtained by replacing a log density with its second order Taylor approximation around the point of maximum. This leads to a Gaussian measure centered at the maximum with a covariance corresponding to the Hessian of the negative log-density (see, e.g., (Bishop, 2006, Section 4.4)). The asymptotic behavior of the parametric Laplace approximation in the small noise or large data limit has been studied extensively in the past (see, e.g., Wong (2001)). The asymptotic approximation of general integrals of the form e λf (x) g(x) dx by Laplace's method is presented in Olver (1974); Wong (2001). Non-asymptotic error bounds for the Laplace approximation can be found in Olver (1968) for the univariate case and in McClure and Wong (1983); Inglot and Majerski (2014) for the multivariate case. Lapiński (2019) studied the Laplace approximation error and its convergence in the limit λ → ∞ in the multivariate case when the function f depends on λ . Coefficients in the asymptotic expansion of the approximated integral are given in Nemes (2013).
Laplace approximation is an important step in establishing the prominent Bernsteinvon Mises (BvM) Theorem that quantifies the convergence of the scaled posterior distribution toward a Gaussian distribution in the large data or small noise limit. Parametric BvM theory is well-understood Van der Vaart (2000); Le Cam (2012). Modern applications with a high dimensional parameter space and limited sample size pose new questions and identify new issues in study of applicability and accuracy of Laplace approximation.
We refer to Lu (2017) for a study a parametric BvM theorem for nonlinear Bayesian inverse problems with an increasing number of parameters. A number of papers discuss the BvM phenomenon for nonlinear inverse problems; see e.g. Nickl (2020); Monard et al. (2019); Giordano and Kekkonen (2020), where the convergence is quantified in a distance that metrizes the weak convergence. Schillings et al. (2020) showed that the Laplace approximation error in Hellinger distance converges to zero in the order of the noise level.
The recent paper Helin and Kretschmann (2022) provides a finite sample error of Laplace approximation for the total variation distance with an explicit dependence on the dimen-sion and of the nonlinearity of the forward mapping. The Laplace approximation is also widely utilized for different purposes in computational Bayesian statistics; see e.g. Rue et al. (2009). 1.1 Motivation: Gaussian approximation of the posterior As one of the main motivation for this study, consider the problem of Bayesian inference for the log-likelihood function L(θ) = L(Y , θ) with data Y , a parameter θ ∈ Ê p and a Gaussian prior π ∼ N (θ 0 , G −2 ) . Here G −2 is a symmetric positive definite matrix in Ê p . Then the posterior density π G (·) of θ given Y is proportional to the product e L(θ) e − G(θ−θ 0 ) 2 /2 : where the sign ∝ means equality up to a normalizing multiplicative constant. Assume that the penalized maximum likelihood estimator (pMLE) θ G is well defined: Clearly θ G is also maximizer of π G (θ) . That is why it is often referred to as maximum a posteriori (MAP) estimator. Let also the log-likelihood function L(θ) be twice differentiable and weakly concave. Define Assuming the latter expression to be positive definite for all considered θ , define also its square root D G (θ) = G (θ) . We use the shortcut D G = D G ( θ G ) . Laplace's approximation means that the posterior distribution π G is close to the Gaussian distribution N ( θ G , D −2 G ) . A closely related Bernstein -von Mises phenomenon claims an approximation of the posterior by N ( θ, D −2 ) , where θ is the MLE and D 2 = = −∇ 2 L(θ * ) is the Fisher information matrix for the true parameter value θ * ; see e.g. Schillings et al. (2020) for a detailed discussion in context on nonlinear inverse problems. The mentioned results provide an efficient tool for Bayesian uncertainty quantification and constructing the elliptic credible sets as level sets of the approximating Gaussian distribution; see Helin and Kretschmann (2022) or Reich and Rozdeba (2020) for applications to drift and diffusion estimation.

Motivation: Gradient free optimization
Suppose that the point of maximum x * of the function ℓ(x) is not known and has to be evaluated numerically. We also assume that the function ℓ(·) is sufficiently smooth and can be efficiently computed at any point x , however, its gradient is not available. Nesterov and Spokoiny (2017) offered a powerful gradient free method based on averaging the exponent e ℓ(x) w.r.t. to some Gaussian distribution for x . Here we aim at reconsidering this idea within the Bayesian optimization framework. Namely, we intend to design a procedure based on Laplace iterations with an updated Gaussian prior at each step. Let π 0 be a starting Gaussian prior π 0 ∼ N (x 0 , G −2 0 ) . The corresponding posterior is defined by normalizing the product e ℓ(x) π 0 (x) . Due to the Laplace approximation result, this posterior is nearly Gaussian with the mean x 1 = argmax ℓ(x) − G 0 x 2 /2 ; see e.g. Theorem 2.1 below. This leads to the general idea of Bayesian optimization: use the posterior mean/variance to update the prior and repeat the Laplace approximation step. One may use standard Monte-Carlo or quasi Monte-Carlo methods for numerical approximation of the posterior mean; see e.g. Schillings et al. (2020) in context of Bayesian inverse problems.

Challenges
In spite of numerous existing results on Laplace approximation, some questions are still open. A small list of relevant issues is given below.
Prior impact in high dimension Modern applications in statistics and optimization force to extend the scope of applicability of the classical results on Laplace approximation by including the cases of a very large dimension even exceeding the sample size. It appears that in the case of a high-dimensional parameter, the prior is not washing out from the posterior. In the contrary to the classical theory, it becomes extremely important. It is of crucial importance to understand well the impact of the prior and its proper choice.
Smoothness conditions Standard smoothness conditions in terms of a uniform bound on the third derivative of f could be too restrictive and hard to check.
The use of posterior mean Another issue is possibility of using a Gaussian approximation of the posterior with inexact parameters. Indeed, the pMLE/MAP θ G is hard to compute, the classical MCMC based Bayesian computations deliver an estimate of the posterior mean and of posterior covariance. A big question is to justify the use of these quantities for uncertainty quantification or Bayesian optimization.
Non-concavity The classical results of Laplace approximation are stated under the condition that the function f (·) is strongly globally concave. In many applications including nonlinear inverse problems Nickl (2020); Schillings et al. (2020); Helin and Kretschmann (2022) or Deep Neuronal Networks such an assumption appears to be unrealistic.

This paper's contributions
This paper aims at reconsidering the classical results on Laplace approximation and to address the above issues. Below the list of the most important achievements in the paper.
Effective dimension and dimension free guarantees We introduce the notion of effective dimension p G of the problem which can be small of moderate even for huge parameter dimension p . The value p G is defined by an interplay between the information delivered by the data and information contained in the prior; see Section 2.1.2 for more details. Further we establish explicit non-asymptotic and dimension free guarantees for the accuracy of Gaussian approximation of the posterior in total variation (TV) distance in terms of effective dimension; see Theorem 2.1. In the case when the nonpenalized log-likelihood function grows linearly with the sample size n , the quality of Laplace approximation is of order p 3 G /n . It can be improved to p 3 G /n if instead of TV-distance, we limit ourselves to the class of centrally symmetric sets. The proofs combine classical variational arguments with sharp bounds for Gaussian quadratic forms.
Conditions require that f is strongly concave and locally smooth with a uniform bound on the third Gateaux derivative of f in a local vicinity of the point of maximum.
Critical dimension The result of Theorem 2.1 helps to address the issue of critical dimension for applicability of Laplace approximation: the relation p 3 G ≪ n between the sample size n and the effective dimension p G is sufficient for our main results. The result on concentration of the posterior only requires p G ≪ n .
Posterior mean in place of MAP The use of the posterior mean x in place of the MAP x * is justified by Theorem 2.10 under the same condition p 3 G ≪ n . The obtained results justify a Bayesian optimization procedure based on iterated Laplace approximation with a sequentially adjusted prior.
Non-concavity Section 3 indicates for the special case of nonlinear inverse problem how the assumption of strong concavity of f (·) can be relaxed using a "warm start" condition which means that a good starting guess is available.
The paper is organized as follows. Section 2 presents general bounds on the error of Laplace approximation and also discusses the use of posterior mean in place of posterior mode as well as the Bayesian optimization procedure based on Laplace iterations. Section 3 specifies the results to the case on nonlinear inverse problem. Proofs and technical tools are collected in the Appendix.

Dimension free bounds for Laplace approximation
Here we present some general results on accuracy of Laplace approximation.

Setup and conditions
Let f (x) be a function in a high-dimensional Euclidean space Ê p such that e f (x) dx = C < ∞ , where the integral sign without limits means the integral over the whole space Ê p . Then f determines a distribution È f with the density C −1 e f (x) . Let x * be a point of maximum: We also assume that f (·) is smooth, more precisely, three or even four time differentiable.
We aim at approximating the measure È f by a Gaussian measure N (x * , D −2 ) . Given a function g(·) , define its expectation w.r.t. È f after centering at x * : (2.1) A Gaussian approximation I D (g) for I(g) reads: The choice of the distance between È f and N (x * , D −2 ) depends on the considered class of functions g . The most strong total variation distance can be obtained as the supremum of |I(g) − I D (g)| over all measurable functions g(·) with |g(u)| ≤ 1 : The results can be substantially improved if only centrally symmetric functions g(·) with g(x) = g(−x) are considered. Obviously Moreover, as x * = argmax x f (x) , it holds ∇f (x * ) = 0 and Implicitly we assume that the negative Hessian D 2 = −∇ 2 f (x * ) is sufficiently large in the sense that the Gaussian measure N (0, D −2 ) concentrates on a local set U . This allows to use a local Taylor expansion for f ( is also strongly concave, then the mass of È f of the complement of U is exponentially small yielding the desirable Laplace approximation.
Motivated by applications to statistical inference, we consider f in a special form for a symmetric p -matrix G 2 ≥ 0 . Here ℓ(·) stands for the log-likelihood function while the quadratic penalty G(x − x 0 ) 2 /2 corresponds to a Gaussian prior N (x 0 , G −2 ) . We also assume that ℓ(·) is concave with D 2 def = −∇ 2 ℓ(x) ≥ 0 . Then clearly In typical asymptotic setups, the log-likelihood function ℓ(x) grows with the sample size or inverse noise variance, while the prior precision matrix G 2 is kept fixed. Decomposition (2.5) is of great importance for obtaining the dimension free results. The main reason is that the quadratic penalty does not affect smoothness properties of the function ℓ(·) but greatly improves the quadratic approximation term. Now we state precise conditions.
More specifically, we assume the following condition.
(C 0 ) There exists an operator G 2 ≤ −∇ 2 f (x * ) in Ê p such that the function If ℓ(·) in decomposition (2.5) is concave then this condition is obviously fulfilled.
More generally, if ℓ(·) in (2.5) is weakly concave, so that ℓ( plays an important role in our conditions and results. Remark 2.1. The condition of strong concavity of f on the whole space Ê p can be too restrictive; see the example of nonlinear inverse problem in Section 3. This condition can be replaced by its local version: there exists a set X 0 such that the Gaussian prior N (x 0 , G −2 ) concentrates on X 0 with a high probability and the maximizer x * G belongs to X 0 . In all the results, the integral over Ê p has to replaced by the integral over X 0 .

Effective dimension
With decomposition (2.6) in mind, we use another notation for D 2 = −∇ 2 f (x * ) : Also we write I G (g) instead of I D (g) in (2.2). The effective dimension p G is given by Of course, p G ≤ p but a proper choice of the penalty G 2 in (2.5) allows to avoid the "curse of dimensionality" issue and ensure a small or moderate effective dimension p G even for p large or infinite. The value p G helps to describe a local vicinity U around x * such that the most of mass of È f concentrates on U ; see Section A.4. Namely, let us fix some ν < 1 , e.g. ν = 2/3 , and some x > 0 ensuring that e −x is our significance level. Define (2.8)

Local smoothness conditions
Let p ≤ ∞ and let f (·) be a three times continuously differentiable function on Ê p . We fix a reference point x and local region around x given by the local set U ⊂ Ê p from (2.8). Also consider the second order Taylor approximation , u ⊗2 and similarly the third order expansion and introduce the remainders The use of the Taylor formula allows to bound It is worth noting that the quadratic penalty and the quadratic term in definition of the values δ k (x, u) cancels, k ≥ 3 . Local smoothness of f (·) or, equivalently, of ℓ(·) , at x will be measured by the value ω(x) : We also denote Our results apply under the condition ω ≪ 1 . Local concentration of the measure È f requires ω ≤ 1/3 ; see Proposition A.11. The main results about Gaussian approximation of È f are valid under a stronger condition ω p G ≤ 2/3 with the effective dimension p G from (2.7).

Error bounds for Laplace approximation
Our first result describes the quality of approximation of the measure È f by the Gaussian measure N (x * , D −2 G ) with mean x and the covariance D −2 G in total variation distance. In all our result, the value x is fixed to ensure that e −x is negligible. First we present the general results which will be specified later under the self-concordance condition.

General bounds
Theorem 2.1. Suppose (C 0 ) . Let also p G be defined by (2.7) and r G and U by (2.8).
Introduce the following conditions.

Corollary 2.3. Under the conditions of Theorem 2.2, it holds for
.

Bounds under self-concordance
Often the smoothness properties of f are described in terms of the third derivative.
Our approach does not require a bounded third or fourth full dimensional derivative, instead we consider directional Gateaux derivatives. Let ℓ(·) be three times continuously differentiable and ∇ 3 ℓ(·) stand for the third derivative. Then the Taylor expansion of the third order implies ω ≤ α/3 with Moreover, in many applications, the function ℓ(·) is of the form −ℓ(x) = nh(x) for a fixed smooth function h(·) satisfying the following condition: for a strongly convex function h(·) and (S 4 ) −ℓ(·) = nh(·) for a strongly convex function h(·) and This is a local version of the so called self-concordance condition; see Nesterov (1988).
Under this condition, we can easily bound the values δ k (x * , u) for k = 3, 4 , ω and α : Finally we state the bounds under (S 3 ) and (S 4 ) .

Critical dimension
Here we briefly discuss the important issue of critical dimension. Theorem 2.1 states concentration of È f under the condition ω ≤ 1/3 and Gaussian approximation under the stronger condition ω p G ≤ 2/3 . Moreover, self-concordance (S 3 ) and Lemma A.12 enables us to bound ω p G /n and ω p G p 3 G /n . Hence, concentration of È f requires the critical dimension condition p G ≪ n , while Gaussian approximation applies under p 3 G ≪ n . We see that there is a gap between these conditions. We guess that in the region n 1/3 p G n , a non-Gaussian approximation of the posterior is possible.
Bochkina and Green (2014) provides examples of non-Gaussian posterior limits for nonregular models.

Kullback-Leibler divergence
Theorem 2.1 through 2.4 quantify the approximation È f ≈ N (x * , D −2 G ) in the total variation distance. Another useful characteristic could be the Kullback-Leibler (KL) divergence between È f and N (x * , D −2 G ) . The KL divergence K (È 1 , È 2 ) = 1 log(dÈ 1 /dÈ 2 ) is asymmetric, K (È 1 , È 2 ) = K (È 2 , È 1 ) with few exceptions like the case of Gaussian measures È 1 and È 2 . Moreover, K (È 1 , È 2 ) can explode if È 2 is not absolutely continuous w.r.t. È 1 . We present two bounds for each ordering. For ease of presentation, we limit ourselves to the case when either (L 3 ) or (S 3 ) meets. (2.18) Now we briefly discuss the value K (È G , È f ) . We already know that È f concentrates on U and can be well approximated by È G on U . However, this does not guarantee a small value of K (È G , È f ) . It can even explode if e.g. È f has a compact support. In fact, the log-density of

Mean and MAP
Here we present the bound on |I(g) − I G (g)| for the case of a linear vector function (2.20) Now we specify the result for the special choice Q = D .
2. An interesting question is whether the result of Theorem 2.7 or Corollary 2.8 applies with Q = D G . This issue is important in connection to inexact Laplace approximation; see the next section. The answer is negative. The problem is related to the last term C e −x in the right hand-side of (2.20). The constant C involves the moments of Q(X − x * ) 2 which explode for Q = D G and p = ∞ .

Inexact approximation and the use of posterior mean
Now we change the setup. Namely, we suppose that the true maximizer x * of the function f is not available, but x is somehow close to the point of maximum x * .
Similarly, the negative Hessian D 2 G (x * ) = −∇ 2 f (x * ) is hard to obtain and we use a proxy H 2 . We already know that È f can be well approximated by N (x * , D −2 ) with . This section addresses the question whether N (x, H −2 ) can be used instead. Here we may greatly benefit from the fact that Theorem 2.1 provides a bound on the total variation distance between È f and N (x * , D −2 Therefore, it suffices to bound the TV-distance between two Gaussian distributions: N (x * , D −2 G ) naturally appears in the Laplace approximation, the second one is a proxy used instead. Pinsker's inequality provides an upper bound: for any two measures P, Q where K (P, Q) is the Kullback-Leibler divergence between P and Q . The KL-divergence between two Gaussians has a closed form: However, Pinsker's inequality is only a general upper bound which applied to any two distributions P and Q . If P and Q are Gaussian, it might be too rough. Particularly the use of tr B 2 is disappointing, this quantity is full dimensional even if each of D −2 G and H −2 has a bounded trace. Also dependence on to the class B el (Ê p ) of elliptic sets A in Ê p of the form A = u ∈ Ê p : Q(u − x) ≤ r for some linear mapping Q : Ê p → Ê q , x ∈ Ê p , and r > 0 . Given two symmetric q -matrices Σ 1 , Σ 2 and a vector a ∈ Ê q , define where λ 1 is the vector of eigenvalues of Σ 1 arranged in the non-increasing order and similarly for λ 2 . Götze et al. (2019) stated the following bound for γ 1 ∼ N (0, Σ 1 ) and Fr /3 , k = 1, 2 . By the Weilandt-Hoffman inequality, λ 1 − λ 2 1 ≤ Σ 1 − Σ 2 1 , see e.g. Markus (1964). Here M 1 = tr |M | = j |λ j (M )| for a symmetric matrix M with eigenvalues λ j (M ) .
Furthermore, for X ∼ È f and γ ∼ N (0, I I p ) , any linear mapping Q : As a special case, we consider the use of the posterior mean Theorem 2.10. Assume the conditions of Theorem 2.7 and Theorem 2.9. Then it holds for any linear mapping Q : The same bound holds with Q = n 1/2 I I p in place of D provided that D 2 ≥ C 0 n I I p for some fixed C 0 > 0 . We may conclude that the use of posterior mean x in place of the posterior mode x * is justified under the same condition on critical dimension p 3 G ≪ n as required for the main result about Gaussian approximation.

Bayesian optimization and iterated Laplace approximation
Suppose that the point of maximum x * of the function ℓ(x) is not known and has to be evaluated numerically. We also assume that the function ℓ(·) is sufficiently smooth and can be efficiently computed at any point x , however, its gradient is not available.
Nesterov and Spokoiny (2017) offered a powerful gradient free method based on averaging the exponent e ℓ(x) w.r.t. to some Gaussian distribution for x . Inspired by the results of previous sections, we propose here a modified version of Nesterov and Spokoiny (2017) based on iterated Laplace approximations. Starting from some π 0 = N (x 0 , G −2 0 ) , we iteratively update the Gaussian prior π k ∼ N (x k , G −2 k ) and try to numerically assess the corresponding posterior obtained by normalization of e f k (x) for f k (x) = ℓ(x) − G k (x − x k ) 2 /2 . This posterior is known to be nearly Gaussian by Theorem 2.1, its mean x k+1 can be efficiently estimated by Monte-Carlo or quasi Monte-Carlo sampling; see e.g. Schillings et al. (2020). We use x k+1 for building the Gaussian prior for the step k + 1 . Moreover, due to our results, the value x k+1 is nearly the maximizer of the corresponding quadratic approximation of f k (x) in the vicinity of x k . One can say that computing the posterior mean x k+1 mimics well a step of the second order Newton-Raphson optimization method but does not require computing the gradient and the inverse Hessian matrix. The prior precision matrix G 2 k can be taken in the form λ −1 k G 2 0 , where λ k replaces the step size, λ k → 0 as k increases. The procedure can be written as follows.
The starting guess is important as well the choice of the step multiplier a . Moreover, one can use the variable multiplier a k and incorporate the posterior covariance for variance reduction schemes. One can make a rigorous analysis of the convergence of this algorithm similarly to Nesterov and Spokoiny (2017).
The proposed procedure is closely related to ensemble Kalman filtering technique; see e.g. Schillings and Stuart (2017), Reich (2022) and references therein.
Algorithm 1 Laplace iterations 1. Start with k = 0 ; Fix x 0 and G 2 0 ; 2. Draw a sample (x (m) ) m≤M from N (x k , G −2 k ) . For each m ≤ M , compute w (m) = e ℓ(x (m) ) and update the sums W k = m w (m) and S k = m w (m) x (m) .
3. Compute the value x k+1 = S k /W k as the next prior mean; update the posterior precision matrix G 2 k+1 = aG 2 k . 4. Increase k by one and repeat from Step 2. 5. Iterate until convergence.

Laplace approximation for non-linear inverse problem
Let m(x) = m i (x), i ≤ n ∈ Ê n be a nonlinear mapping (operator) of the source signal x ∈ Ê p to the target space Ê n . We consider the problem of inverting the relation z = m(x) : given an image vector z ∈ Ê n , recover the corresponding source x ∈ Ê p . This leads to the nonlinear least square problem of maximizing the function ℓ(·) of the form Given a Gaussian prior N (x 0 , G −2 ) , we consider Laplace's approximation for the penalized function Laplace's approximation requires weak concavity of ℓ(x) ; see However, such a condition seems can be hard to ensure for all x ∈ Ê p unless m(·) is linear. Instead we restrict ourselves to some subset X 0 ⊂ Ê p which can be viewed as a concentration set of the Gaussian measure N (x 0 , G −2 ) . As the prior mass of the complement of X 0 is exponentially small, when sampling from the prior N (x 0 , G −2 ) , one would need exponentially many samples to hit any set A outside of X 0 . Therefore, we restrict the prior to the set X 0 by skipping all draws of x with x ∈ X 0 . In what follows we also assume that x 0 is a reasonable guess ensuring that the target x * G belongs to X 0 . Then the result about Laplace approximation applies even after restricting the parameter space to X 0 . This local set X 0 is defined as elliptic concentration set for the Gaussian prior N (x 0 , G −2 ) in the form where Q is a linear operator in Ê p and r 0 is fixed to ensure that the most of prior mass is within X 0 . Theorem B.4 suggests X 0 in form of a ball around x 0 with the radius of order √ tr G −2 . A proper choice of G 2 ≥ I I p has to ensure a small value of tr(G −2 ) . Note that the operator Q can be scaled by any factor. Then the radius r 0 in (3.5) will be scaled correspondingly.
Introduce for any x ∈ X 0 a p -symmetric matrix This is the first term in expansion (3.3) corresponding to a linear approximation of m(·) at x . Injectivity of m(·) means thatD 2 (x) is positive definite and well conditioned. If n −1D2 (x 0 ) ≤ G 2 then the choice Q =D(x 0 ) in (3.4) is our alternative to Q = I I p .
Below we assume the following regularity condition.
(∇m) Let the local set X 0 be defined by (3.4) and (3.5) with Q =D 0 def =D(x 0 ) . For some fixed C 2 , C n , and any x ∈ X 0 , u ∈ Ê p , it holds

7)
and with some C 0 ≥ 1D Condition (∇m) only requires some local regularity of the m i (x) 's and injectivity of the gradient mapping ∇m(x) within X 0 . The shape of X 0 is defined by the prior covariance G −2 which has to be selected to ensure that the set X 0 is a local set. More precisely, we connect the size of X 0 and local regularity of m(·) by the following condition.
(r 0 ) With X 0 defined by (3.4) and (3.5) and C n , C 0 , and C 2 from (∇m) , it holds To get a feeling of this condition, consider the case withD 2 = nI I p . Then (3.5) yields r 2 0 ≍ n tr G −2 and (3.9) requires tr G −2 ≪ 1 , that is, the prior N (x 0 , G 2 ) is supported to a small vicinity of x 0 . Now we introduce the key condition of "warm start" which informally means a good choice of the starting point x 0 ensuring that the local vicinity X 0 of x 0 contains a point Later we show that (∇m) , (r 0 ) , and (x 0 ) ensure for some δ ≤ 3/4 , and hence, ℓ(x) is strongly concave for x ∈ X 0 . Now we consider the point x * G from (3.2) and its local vicinity. It is important to secure that x * G is within X 0 . This appears to be automatically fulfilled if the true solution x * = m −1 (z) is within X 0 ; see Lemma A.14. Now we are back to the setup with a strongly concave function ℓ(x) for x ∈ X 0 . Similarly to the general case, define The effective dimension is given by p G = tr(D 2 D −2 G ) ; see (2.7). The local vicinity U G of x * G can be defined by the rule (2.8). Alternatively one can useD 2 =D 2 (x * G ) instead of D 2 : The radius r has to be adjusted to ensure that the such defined vicinity of x * is not smaller than U from (2.8): U G ⊆ U . Relation (3.10) suggests using r = 2r G .
For applying the general results of Theorems 2.1 through 2.4, the function ℓ(x) = − m(x)−z 2 /2 has to be sufficiently smooth. We now present some sufficient conditions in terms of the functions m i (x) . These conditions extend (∇m) .
(∇ 4 m) For some C 4 , it holds uniformly over x ∈ U and u ∈ Ê p We are now prepared to state the main result about Laplace's approximation which is a straightforward application of Theorem 2.1 to the measure È f restricted to X 0 .

A Tools and proofs
Here we collect the proofs of the main results and some useful technical statements about the error of Laplace approximation. Below we write x instead of x * . Note that after passing to representation (2.3), many results below apply to any x , not necessarily for x = x * . We only use D 2 G = −∇ 2 f (x) and ω instead of ω(x) . Everywhere we assume the local set U to be fixed by (2.8). We separately study the integrals over U and over its complement. The local error of approximation is measured by As a special case with g(u) ≡ 1 we obtain an approximation of the denominator in (2.3).
In addition, we have to bound the tail integrals The analysis will be split into several steps.

A.1 Overall error of Laplace approximation
First we show how to seam together the error ♦ of local approximation and the bounds for the tail integrals ρ and ρ G ; see (A.2).
Proposition A.1. Suppose that for a function g(u) with g(u) ∈ [0, 1] and some ♦, ♦ g Then with ρ and ρ G from (A.2) Similarly for g(u) ≥ 0 Putting together all these bounds yields (A.5).
The next corollary is straightforward.

A.2 Lower and upper Gaussian measures
This section introduces the lower and upper Gaussian measure which locally sandwich the measure È f using the decomposition from condition (C 0 ) . Denote −∇ 2 f (x) = D 2 G . Definition (2.9) enables us to bound with ω = ω(x) yielding two Gaussian measures which bounds È f locally from above and from below.
The next technical result provides sufficient conditions for their contiguity.

A.3 Local approximation
This section presents the bounds on the error ♦ of local approximation (A.1). The first result only uses ω p G < 1 while the second one also assumes (L 3 ) . These two results allow to bound the total variation distance between È f and N (x, D −2 G ) . The third result is finer and is based on (L 4 ) . It allows to improve the error term of Gaussian approximation over the class of centrally symmetric sets. We also present some extensions for the moments of È f . Proposition A.5. Let ω = ω(x) from (2.9) and p G from (2.7) satisfy Then for any function g(u) with |g(u)| ≤ 1 Moreover, under (L 3 ) , the bound applies with ♦ = ♦ 3 ; see (2.14).

In view of
for a standard normal γ , we derive yielding the second statement.
The result can be extended to the case with a homogeneous function g(u) , e.g.
Proposition A.6. Suppose the conditions of Proposition A.5 and (L 3 ) . Then for m ≥ 1 and any m -homogeneous function g(·) with g(tu) = t m g(u) Proof. Similarly to the proof of Proposition A.5, under (2.12) for R(t) from (A.14) This yields (A.17).
Important special cases correspond to m = 1 .
Proposition A.7. Suppose the conditions of Proposition A.6. Then with γ G ∼ N (0, D −2 G ) , it holds for any linear mapping Q : Ê p → Ê q and any vector a ∈ Ê q Now we state a sharper result based on (L 4 ) .
Proposition A.8. Suppose the conditions of Proposition A.5 and (L 4 ) . Then for any function g(u) with |g(u)| ≤ 1 and g(u) = g(−u) where for a Gaussian element γ G ∼ N (0, D −2 If the function g(·) is not bounded by one but it is symmetric and 2m -homogeneous, i.e. g(tu) = t 2m g(u) , then (A.19) still applies with Proof. Below we write f (3) instead of ∇ 3 f (x) and δ k (u) in place of δ k (x, u) , k = 3, 4 .

A.4 Tail integrals
In this section we also write x in place of x * . Below we evaluate ρ from (A.2) which bounds the integral of e f (x;u) over the complement of the local set U of a special form U = u : Du ≤ ν −1 r G for D from (C 0 ) . Our results help to understand how the radius r G should be fixed to ensure ρ sufficiently small. The main tools for the analysis are deviation probability bounds for Gaussian quadratic forms; see Section B.1.

This yields by r
Integration by parts allows to represent the last integral as This completes the proof of the result (A.22). The second statement (A.23) is about Gaussian probability È T γ ≥ r G for a standard normal element γ , and we derive and (A.23) follows.
The next result extends (A.22).
Proposition A.10. Assume the conditions of Proposition A.9 with for some m ≥ 0 . Then (A.22) can be extended to Proof. The case m > 0 can be proved similarly to m = 0 using m log z ≤ mz .

A.5 Local concentration
Here we show that the measure È f well concentrates on the local set U from (2.8).
Again we fix x = x * .
Proposition A.11. Assume ω ≤ 1/3 . Then Moreover,  A.7 Proof of Theorem 2.4 We make use of the following lemma.
A.8 Proof of Theorem 2.5 and Theorem 2.6 Define C G def = log e − D G u 2 /2 du − log e f (x * ;u) du .
For u = x − x * , it holds as in (2.3) Similarly to (A.15) and (A.16), we can bound On the complement u ∈ U c , we use that The last inequality is based on concavity of ℓ(·) and local approximation ℓ(x * ; u) ≈ − Du 2 /2 within U yielding ℓ(x * ; u) < 0 for u ∈ U c . By Proposition A.10 with We conclude that Similarly This yields e C G ≤ 1 + ♦ 3 + 2e −x and C G ≤ ♦ 3 + 2e −x . Putting all together results in A.10 Proof of Theorem 3.1 The basic idea of the proof is to reduce the statements of the theorem to the case of Theorem 2.1. Conceptually, the most important step of the proof is to show that the function f (x) is strongly concave on the subset X 0 which is a concentration set of Lemma A.13. Suppose (∇m) , (r 0 ) , and (x 0 ) . Then for all x ∈ X 0 Given x ∈ X 0 , denote u = x − x * . The definition of X 0 implies D 0 u ≤ 2r 0 . The first order Taylor expansion yields Hence, by (3.7) and (3.8) The next step is to show that under (x 0 ) , the point x * G is also within X 0 .
Bound (A.31) means that (2.16) is fulfilled with c 3 = 4C 3/2 G C 3 C n . This implies (3.12) similarly to Theorem 2.4. Moreover, (A.32) implies (S 4 ) with c 4 = 8C 2 G C 4 C n , and the second statement of the theorem follows as well.
Now we consider the case of a non-centered quadratic form Bγ, γ /2 + A, γ for a fixed vector A .

B.2 Deviation bounds for Gaussian quadratic forms
The next result explains the concentration effect of Bξ, ξ for a centered Gaussian vector ξ ∼ N (0, V 2 ) and a symmetric trace operator B in Ê p , p ≤ ∞ . We use a version from Laurent and Massart (2000). For completeness, we present a simple proof of the upper bound.


Intriguing Properties of Compression on Multilingual Models

Multilingual models are often particularly dependent on scaling to generalize to a growing number of languages. Compression techniques are widely relied upon to reconcile the growth in model size with real world resource constraints, but compression can have a disparate effect on model performance for low-resource languages. It is thus crucial to understand the trade-offs between scale, multilingualism, and compression. In this work, we propose an experimental framework to characterize the impact of sparsifying multilingual pre-trained language models during fine-tuning. Applying this framework to mBERT named entity recognition models across 40 languages, we find that compression confers several intriguing and previously unknown generalization properties. In contrast to prior findings, we find that compression may improve model robustness over dense models. We additionally observe that under certain sparsification regimes compression may aid, rather than disproportionately impact the performance of low-resource languages.

Introduction
Scaling language models benefits multilingual settings, since it is difficult to maintain performance across a growing number of languages at a constant model size, a property also called the "curse of multilinguality" (Conneau and Lample, 2019;Conneau et al., 2020;Artetxe and Schwenk, 2019). However, the extent of growth in language model (LM) size (Radford et al., 2019;Brown et al., 2020;Zhang et al., 2022;Chowdhery et al., 2022) has made deployment to resource-constrained environments much more challenging (Warden and Situnayake, 2019;Samala et al., 2018;Treviso et al., 2022). To benefit from the performance gains conferred by scale, efficiency techniques that reduce model size while maintaining comparable aggregate performance are widely used, such as quantization (Shen et al., 2020), compression (Michel et al., 2019Lagunas et al., 2021) and distillation (Tsai et al., 2019;Sanh et al., 2019;Pu et al., 2021).
While most compression techniques have minimal impact on aggregate performance numbers (Gale et al., 2019;Li et al., 2020;Hou et al., 2020;Chen et al., 2021;Bai et al., 2020;ab Tessera et al., 2021), the impact on individual sub-populations in the data, such as low-resource languages, can be far more severe (Hooker et al., 2019;Hooker et al., 2020;Ahia et al., 2021). Disparities in resource availability become more apparent at larger scale, both in terms of data and deployment resource availability. This makes compression all the more necessary, but also motivates a thorough consideration of the subsequent impact of compression on generalization.
In this work, we develop an experimental framework to investigate the impact of compression during fine-tuning of pre-trained multilingual models which we apply to Named Entity Recognition (NER) across 40 languages of the WikiAnn benchmark (Pan et al., 2017). We study the impact of compression on groups of languages across multiple dimensions-resourcedness, script, and language family-and evaluate the sensitivity of models to input perturbations along these groupings.
This leads us to discover the following intriguing properties: (1) Lower-performing languages disproportionately suffer under extreme levels of sparsity, as pruning amplifies disparities. However, low-resource languages present an intriguing flipflop moment, where their performance may benefit from medium regimes of sparsity. (2) We find that dense models overfit to typical test cases, achieving a close-to-0 F1 score on slightly perturbed inputs, while compression can recover close to the original test performance. Our results stand in contrast to previous work that find that sparsity erodes robust-ness, suggesting more work is needed to understand the dynamics between compression and robustness.
(3) The choice to prune model embeddings can completely negate the two benefits described in the previous observations, showing the importance of comparing the two cases in future analyses.

Related Work
The "curse-of-multilinguality" creates a tradeoff between number of languages and size of a model (Conneau et al., 2020). However, training smaller models supporting fewer languages may not always be feasible (Abdaoui et al., 2020). Compressing large models has been shown to combat the curse, either by compressing the pretrained model (Tsai et al., 2019;Sanh et al., 2019), or by compressing during finetuning, as in our case. While many studies investigate the impact of pruning on aggregate metrics in monolingual pretrained LMs (Sanh et al., 2020;Goyal et al., 2020;Gordon et al., 2020;Budhraja et al., 2020;Sajjad et al., 2020;Lagunas et al., 2021;Xu et al., 2021;Du et al., 2021a;Ganesh et al., 2021), fewer works focus on multilingual settings (Mukherjee and Hassan Awadallah, 2020;Ansell et al., 2022). Yet, prior analyses find a disparate effect of removing attention heads or model layers on languages and language families distant from the training data in NER (Ma et al., 2021;Budhraja et al., 2021), demonstrating the importance of looking into subpopulations as we do in this study. Concurrent work (Mohammadshahi et al., 2022b) explores compression on multilingual machine translation models, and find that while under-represented languages suffer extreme performance drop, some mediumresource language benefit from compression. However, unlike our proposed framework, they perform this study without additional fine-tuning after compression.
Studies that compare the robustness of compressed and dense models further find that compression may lead to erosion of performance on "challenging" samples and poor generalization (Ahia et al., 2021;Du et al., 2021a;Xu et al., 2021), a finding that we expand on and connect to language resourcedness. The technique we use to study robustness expands on studies that perturb training (Yaseen and Langer, 2021;Dai and Adel, 2020) or evaluation data (Dhole et al., 2021) in NER by introducing perturbations specific to languages, language families, and scripts.

Methodology
Data We conduct our experiments on WikiAnn (Pan et al., 2017), a multilingual NER dataset. WikiAnn was sourced from Wikipedia articles and automatically annotated with LOC (location), PER (person), and ORG (organisation) labels in the IOB2 format (Ramshaw and Marcus, 1995). It is considered a "silver standard" due to its automatic entity labels and noise (Lignos et al., 2022), but with its 176 languages it covers the most languages of any NER dataset. We focus our experiments on the 40 languages from the XTREME benchmark (Hu et al., 2020), with train-test splits defined by Rahimi et al. (2019). These training sets were built with stratified sampling to create a balance across entity types (Lignos et al., 2022), and are thus a subset of the total available data from the original WikiAnn. Table 1 lists language codes in ISO 639-1 and their available training data for fine-tuning.
Perturbations We test the robustness of compressed models by perturbing named entities in the test set. Previous work (Du et al., 2021a) show that sparse pretrained language models are less robust than their dense equivalents when evaluated on adversarial test sets, even when they perform similarly on in-distribution test sets. We adopt a data perturbation technique from Dai and Adel (2020) called entity mention replacement; an entity is randomly swapped with another entity of the same type (example sentences shown in App. D). We first perturb entities within same language for all the languages in our dataset (in-language). Secondly, we propose a new benchmark appropriate for testing the cross-lingual robustness of multilingual models on our downstream task. We perturb entities across different languages that share common linguistic properties. In particular, we group languages by family and script and perturb entities across languages within the same group (in-script, in-family).
Model We use the cased multilingual BERT (mBERT) (Devlin et al., 2019) for all our experiments because it is one of the most widely used and studied multilingual LMs (e.g., Pires et al., 2019;Rönnqvist et al., 2019;Wu and Dredze, 2019;Wang et al., 2020;Chi et al., 2020 ,he,hu,id,it,ja,ko,ms nl,pt,ru,th,tr,ur,vi,zh  trained on Wikipedia data from 104 languages, has approximately 177 million parameters, and a vocabulary size of around 120,000. We finetune mBERT by appending a linear classification layer to the model and updating all its parameters. Full hyperparameters are listed in App. A.1. We evaluate the impact of sparsity in two settings: 1) In the monolingual setting we fine-tune on individual languages. For this setting, we select 10 languages with different available data size, language family and scripts (see Table 1). 2) In the multilingual setting we jointly fine-tune on all languages. We train all models with three random seeds, and evaluate F1 using seqeval (Nakayama, 2018) on the individual languages' evaluation data. We report mean results across runs after computing the micro-average F1 scores across entity classes.
Pruning We induce sparsity by applying Iterative Magnitude Pruning (IMP) (Han et al., 2015(Han et al., , 2016 during fine-tuning. IMP iteratively removes weights that are below a certain threshold until a desired target sparsity is reached. IMP is widely used and competitive with far more compute intenparameters that replicate the prior results.  sive approaches (Gale et al., 2019;Gordon et al., 2020;Du et al., 2021b;Ganesh et al., 2021), while allowing us to sparsify to an exact level. We compare two pruning strategies: 1) partial where we prune all dense layers except for embedding layers, 2) incl. embeddings where we prune all dense weights including embedding layers. Embeddings make up more than half (91M) of the 177M parameters in mBERT, while dense weights make up the rest. Hence, pruning embeddings allows us to significantly reduce the number of mBERT parameters. We consider five sparsity levels: 50%, 70%, 80%, 90%, 95% and 98%, corresponding to the percentage of weights pruned (hyperparameters in App. A.2. Preliminary experiments were conducted with lower sparsity levels (10%-40%) and yielded similar findings to those at moderate sparsity levels (50%-70%), motivating the sparsity intervals chosen. The chosen sparsity levels also align with general best practice in sparsity evaluation as presented in previous works. Moderate to high sparsity levels (50%+) are necessary for efficiency gains in the real-world and are usually studied in literature ( Sparsity (Pruning with Embeddings)

Type of Perturbation
In-Language In-Script In-Family Figure 3: Regular vs Perturbed: We show the aggregated results across all languages after perturbation at different sparsity levels. Without pruning, the model performs poorly, which is overcome by partial pruning, but not pruning with embeddings. The relative performance drop is consistent across all pruning levels above 0.
to benefit more from crosslingual transfer. 2 We find that this finding holds under sparsity -multilingual models achieve higher F1 than monolingual models not only in the dense setting, but across all sparsity levels, as shown in Figure 1. 3 At high sparsity levels, the loss in quality that is generally incurred is considerably lower for multilingual models. This suggests that when high levels of compression are necessary (e.g. for inference efficiency needs), multilingual training should be preferred to monolingual training, as it could help offset some of the erosion in the performance caused by the compression. Thus, we conclude that the benefits of cross-lingual transfer are not inhibited by pruning, and perhaps are even more pronounced at a lower capacity (Dufter and Schütze, 2020) for certain languages.

Impact of pruning across languages
Figure 2 displays the relative differences in F1 score between dense and sparse models across languages, grouped according to fine-tuning size. 4 At moderate sparsity levels (50%-70%), partial pruning surprisingly improves over the dense models, in particular those with less fine-tuning data. The majority of languages (26 out of 40) benefit from moderate pruning and yield slightly higher F1 with pruning than without. All three datasets with only 100 fine-tuning examples (yo, my, jv) benefit. This suggests that moderate pruning may benefit low-resource datasets when introduced during a finetuning regime. However, at high sparsity levels (70%-98%), the findings reverse. Those languages that have a lower frequency of representation in the finetuning dataset incur the highest absolute and relative loss in quality. We can observe the same trend when grouping languages according to their family or script, respectively (see Fig. 4 and 5 in App. B). The groups that start with the lowest average performance under the dense model, also suffer the most under extreme sparsity.
In conclusion, moderate pruning levels should be explored for low-resource languages since they may benefit such languages. This is especially important since models for low-resourced languages are often deployed in resource-constrained environments § (Ahia and Ogueji, 2020;Nekoto et al., 2020;Ahia et al., 2021;Ogueji et al., 2021;Mohammadshahi et al., 2022a). Also, since high sparsity levels reinforce existing disparities (as measured by model performance and data availability) between languages and language groups, it is imperative that practitioners pay attention to possible disparities when sparsifying models. Figure 3 shows the relative performance on the perturbed sets as a fraction of the corresponding unperturbed performance. Across all perturbation types, the dense model performs poorly, indicating that the model may have overfit to typical entities and the semantic context that appear in the training corpora. 5 Surprisingly, partial pruning at any level (2021) find that magnitude pruning improves model robustness to out-of-distribution shifts in machine translation. Despite the contradictions, our work represents an important step in understanding the impact of pruning on robustness, especially since we are one of the firsts to explore it multilingually. Interestingly, our findings are consistent across all perturbation types as their scope increases from languages (in-language) to scripts (in-script) and families (in-family). This suggests that sparsity can be explored as an avenue to improve robustness as has been explored in previous works (Xu et al., 2021;Ahia et al., 2021).

How does pruning impact robustness?
However, pruning the embeddings makes a crucial difference for the perturbed test cases. While pruning the embeddings does not matter for regular test set (see Figure 2), we observe the same severe drop in performance on the perturbed test-set as for the dense model. This suggests that including model embeddings when pruning sharply erodes performance on out-of-distribution rare artefacts, prompting a closer look into what is pruned in the embedding space and the potential impact of sparsifying different parts of a model.

Conclusion
This work investigates the effects of compression on multilingual pre-trained language models during fine-tuning. Our analysis revealed several intriguing properties of pruning that should inform future work in this direction: (1) Pruning dense layers up to ∼ 70% may improve quality for lowfrequency examples in the data and enhance model robustness.
(2) The decision to prune embeddings may have critical impact on model robustness to out-of-distribution performance. (3) While lowperforming languages benefit from moderate pruning, they are disproportionately harmed when pruning more aggressively. Based on these intriguing properties, we also make several recommendations to machine learning practitioners.
performance of dense models on the perturbed test sets.

Limitations
We detail the following potential limitations of our work: Noisy dataset: Lignos et al. (2022) shed light on several quality issues of the WikiAnn dataset that we are treating as a gold standard. Our results might thus not adequately reflect NER performance that can be achieved with cleaner and humanannotated datasets, such as the MasakhaNER (Adelani et al., 2021) or SADiLaR (Eiselen, 2016). Since the perturbations are based on the WikiAnn labels, we might be amplifying the existing label noise for the perturbed test sets and as a result underestimate model quality on clean perturbed examples. We try to combat the randomness by averaging results across three separate runs, but any issues intrinsic to WikiAnn will likely impact all three.
Other Multilingual Models and Downstream tasks: Multilingual pre-trained models such as XLM-R (Conneau et al., 2020) might yield a better performance or show slightly different trends across languages (Adelani et al., 2021). Other downstream tasks, especially generation tasks, might tolerate different levels of sparsity, and also show different crosslingual transfer capabilities (Wu and Dredze, 2019; Hu et al., 2020). However, since fine-grained prior results on the same WikiAnn splits were not available to us, we restricted the analysis to mBERT where we could verify that we can replicate the results reported by XTREME.
Evaluation metrics: We use F1 as the sole evaluation metric and trust it to reflect quality adequately across languages. Human evaluation and the use of qualitative evaluation metrics might reflect the quality for individual languages better.

Unknown factors influencing performance:
The absolute performance for a given language can be influenced by many factors including size, family and script, relatedness to other languages, and the inherent difficulty of the NER task and the evaluation examples, as studied in related works (e.g., Pires et al., 2019;Wu and Dredze, 2020;Shaffer, 2021;Adelani et al., 2021;Muller et al., 2021;Deshpande et al., 2021). As a result, it is impossible to identify the exact cause for all our observations and we have to partially rely on correlational analyses.  Figure 9 shows the correlation between percentage entity overlap and F1 on dense multilingual models.

C Full Results
We present the results for individual languages on both the regular and perturbed test sets obtained via multilingual finetuning in tables 2, 3, 4 and 5 We present the results for individual languages on both the regular and perturbed test sets obtained via monolingual finetuning in tables 6, 7, 8 and 9

D Examples of Perturbed Test Sentences
We present examples of perturbed test sentences in the in-language setting for English (table 10) Table 9: F1 scores for monolingual fine-tuning on the perturbed data for various levels of sparsity with pruning embedding layers.

Example English test sentences
Original Much construction was undertaken during this period , such as the building of Palermo Cathedral . Perturbed Much construction was undertaken during this period , such as the building of Knott 's Soak City .

Original
It is found in Peru . Perturbed It is found in Carbon Cliff , Illinois .

Original
It flows from Ägerisee through Lake Zug into the Reuss . Perturbed It flows from New Orleans through Humboldt County , Nevada into the Crow Agency , Montana .

Original
The album 's lead single " Better Believe It " featuring Young Jeezy and Webbie , was released on July 14 , 2009 . Perturbed The album 's lead single " Better Believe It " featuring W. S. Merwin and Empress Maria Theresa , was released on July 14 , 2009 . ' " ' " jé Fáráò ni Égíptì Ayéijóun . Perturbed ' " ' " j Yousaf Raza Gillani ni Nàìjíríà . Table 11: Example of test sentences for Yoruba language using the entity mention replacement (Dai and Adel, 2020) technique where an entity is randomly swapped with another entity of the same type.


A Fuzzy Crew Rostering Model Based on Crew Preferences and Seniorities considering Training Courses: A Robust Optimization Approach

Crew scheduling problem is divided into crew pairing problem (CPP) and crew rostering problem (CRP). In this paper, a rostering model is presented to assign crew to pairings in such a way that total weighted preference is maximized. Crew members declare which parings they wish to be assigned and which ones are undesirable for them. A score is calculated in the objective function if a crew member is assigned to his/her preferred pairing, and a penalty is considered if he/she is assigned to an undesirable pairing. Moreover, crew seniorities are considered in calculating total preference. In addition, the model considers standard rules and regulations as well as crew attendance at the required training courses. The model is formulated in such a way that inconsistent crew members are not assigned to a flight. Due to the uncertainty in determining of the seniority weight, this parameter is considered as fuzzy. At the end, the robust counterpart of the nominal model is developed due to the uncertainty of time away from the base (TAFB). In this research, the issue of inconsistent crew in rostering problem is considered for the first time. Moreover, a new scoring mechanism is introduced to calculate desirable and undesirable assignments in the objective function. The proposed CRP is solved using the genetic algorithm (GA), and its performance is verified in comparison with GAMS in some test problems. On average, the optimality gap in GA is only 0.5 percent. Finally, the proposed model is examined with real-world data from Air India Airline. In comparison with the previous research studies, the suggested model (scoring mechanism) reduced the number of undesirable rosters by 61.59%.

Introduction
Using the mathematical models in airline operations is an interest field for researchers. Because of the complexity of optimization problems in the airlines, they are usually divided into several smaller problems: flight scheduling (FS), fleet assignment (FA), aircraft maintenance routing (AMR), and crew scheduling (CS) [1,2]. Flight scheduling is the starting point of airline operations [3]. Flight scheduling output is utilized as fleet assignment input. e purpose of fleet assignment is to assign fleet types to flights scheduled in the previous step. Following FA problem, individual aircraft routing is planned so that maintenance requirements are satisfied [4].
One of the major topics in airlines is crew scheduling problem (CSP). Due to its complexity, CSP is separated into two phases: crew pairing problem (CPP) and crew rostering problem (CRP) [5]. Pairings are generated in CPP based on carrier rules and regulations [6]. A pairing includes a set of sequential duty days separated by layovers. e purpose of CPP is to cover all flights in such a way that crew cost is minimized. e pairings created in the CPP are assigned to crew in CRP, and each (co)pilot's roster is determined [7]. Although the CPP focuses on cost minimization, the aim of CRP is to maximize crew satisfaction by assigning balanced workload and meeting crew preferences.
ere are three different approaches in CRP: bidline, personalized with strict seniority, and personalized with global objective [8]. In bidline approach, which was common in North American airlines, anonymous monthly schedules, called rosters, are generated. en, crew members bid for the rosters according to their preference based on seniority. In a personalized approach with strict seniority, which is an emerging approach in North America, crew preferences are met sequentially in decreasing order of seniority. In personalized approach with global objective, crew assignment is done in such a way that total satisfaction is maximized without considering seniority [9].
Crew rostering problem is one of key the phases of airline scheduling. If crew members are assigned to their favorite rosters, a high degree of satisfactions is resulted. Using mathematical models in rostering problem leads to crew satisfaction. In this paper, crew members are assigned to rosters in such a way that maximum satisfaction is reached. To calculate the satisfaction, a new scoring system is used. On the other hand, the issue of travel safety is one of the most important issues for airlines. Assignment of inconsistent crew to a single flight threatens flight safety considerably. As a main motivation, we develop the proposed model to prevent assignment of inconsonant crew to a single flight.
A new model for CRP is formulated in this paper. e objective function maximizes total weighted preference. Instead of using simple bidline or personalized approach, a new mechanism for calculating preference is used in this paper. At first, each crew member declares if he/she prefers to be assigned to a given pairing or not. If they are assigned to one of their preferred pairings, one positive score is calculated in the objective function. Similarly, a negative score (penalty) is considered if a crew member is assigned to an undesirable pairing. e objective function maximizes crew preferences considering seniority. On the other hand, the model considers a crucial qualitative factor in crew assignment to flights. As some of the crew members may be inconstant and have bad work relations to each other, the model prevents to assign them to a single flight. Assignment of the inconsistent crew to a flight may lead to irrecoverable damage. Moreover, crew attendance at the required training courses and some of the airline regulations such as maximum time away from the base (TAFB) and maximum flying hours are considered in assigning crew to rosters. As TAFB is an uncertain parameter in real-world cases, robust counterpart of the proposed model is formulated using the approach of Bertsimas and Sim [10]. A genetic algorithm is used to solve the model for an Indian airline (Air India). GA efficiency in solving the proposed model is examined by an exact solution obtained by CPLEX in GAMS for a variety of small and medium-size examples. is paper benefits from two main contributions. e first one is using a scoring mechanism instead of bidding for preferred pairings. Against previous research studies, both willingness and unwillingness of crew in assigning to each pairing are considered. e second one is considering the issue of inconsistent crew. is important qualitative factor should be considered to ensure a safe travel for passengers. e organization of the paper is as follows. e literature on the crew rostering problem is discussed in Section 2. In Section 3, problem statement and mathematical model are described in detail. Section 4 discusses the solution approach, application of the CRP model in a real-world case study, and the obtained results. Finally, Section 5 concludes with findings and ideas for future studies.

Literature Review
In this section, a comprehensive review on the CRP is presented. Different research studies have been classified based on if their contributions are in solution approach (heuristic, metaheuristic, and exact) or development of the mathematical model. Although many studies have concentrated on developing solution approaches, only a few ones have contribution in mathematical formulation.

Development of the Solution Approach.
Boubaker et al. [11] proposed two heuristic algorithms to solve the CRP. e first approach was the standard branch-and-price algorithm, and the second one was a combination of Elhallaoui et al. [12] with the branch-and-price. Computational findings indicated that the combined method for large-scale problems performed better than standard ones. Maenhout and Vanhoucke [13] considered the use of the additional crew and developed a scatter search algorithm to solve the CSP. e result of comparing the proposed approach with the branch-and-price algorithm and variable neighborhood search demonstrated that this technique could result in highquality solutions. Saddoune et al. [14] discussed the weaknesses of the traditional approaches, solved the CRP using the dynamic column generation method, and implemented it in an American airline. An integrated approach for aircraft and crew schedule recovery problems was proposed by Zhang et al. [15]. To solve the real problem, they used a twolevel metaheuristic algorithm and showed that it could lead to high-quality solutions. Zhang et al. [16] used the Tabu search algorithm to solve the CRP and used the concept of search memory, which accelerated the convergence. By testing the samples, the findings indicated that their method achieves near to optimal solutions in a reasonable time. Doi et al. [17] addressed obtaining a fair roster for the crew by applying a two-level decomposition approach. Zhang et al. [18] developed three multiobjective evolutionary algorithms to search for Pareto solutions. e experimental results on a Chinese airline proved that the proposed technique was capable to improve the rosters' quality by balancing the crew workload. Zhang et al. [19] presented an implicit mathematical model using bidline approach and developed a new variable neighborhood search to solve it. e results indicated that the proposed approach could provide high-quality solutions for two different scenarios of the crew workload. Chutima and Arayikanon [20] formulated a CSP in a low-cost airline and solved it using a multiobjective evolutionary optimization algorithm (MOEA). Comparing the proposed algorithm results with the honeybee mating optimization (HBMO) indicated that MOEA was superior in terms of convergence criteria. e improved discrete particle swarm optimization (IDPSO) method was proposed by Zheng [21] to solve CSP. e findings showed that IDPSO's performance was very efficient. Zhou et al. [22] improved an ant colony algorithm to solve the CRP in a biobjective model optimizing fairness and satisfaction. ey proposed a hybrid complementary heuristic strategy and a local search strategy to solve the model adequately. ey showed that the novel method led to higher-quality solutions than the greedy algorithm and other prominent multiobjective optimization techniques. Chutima and Krisanaphan [23] considered crew cost minimization as well as balancing the workload of pilots. ey solved the proposed model with the adaptive nondominated sorting differential evolution algorithm III (ANSDE-III) and compared the results with the NSGA-III and the multiobjective differential evolutionary approach (MODE).
e results showed that ANSDE-III performed better than other algorithms. [24] formulated pilot roster recovery problem with the aim of minimizing the maximum flight time of each pilot after a disruption and reduced exchange costs. ey solved the model with the nondominated sorting genetic algorithm II (NSGA-II) to find Pareto solutions and applied the constraint-loosening mechanism (CLM). e results indicated that the proposed method led to high-quality solutions to the CR problem. Crew changeability in case of delay and disruption was offered by Ilagan and Sy [25]. According to the results, although the proposed model had increased costs, the delay had been significantly reduced. Kasirzadeh et al. [26] formulated the personalized CR problem with undetermined pairings and fixed pilots. e model aimed to optimize crews' preferences and their cost. ey used the approach of Saddoune et al. [14] and obtained an acceptable level of crew satisfaction. Armas et al. [27] proposed the CRP with the aim of fair distribution of workload among crew members, and they applied regulations that had not been considered in previous studies. To solve this model, they used the multistart randomized heuristic algorithm and showed that this algorithm has the ability to solve the problem with high quality. An integrated model for CSP was presented by Zeighami and Soumis [8].

Development of CR Models. Chen and Chou
ey considered pilots' vacation requests in the CP to generate better pairings.
e authors formulated the model by minimizing pairing costs and maximizing the desired number of vacations.
e experiment results showed that combining Benders' decomposition algorithm and column generation could significantly save CPU time. Quesnel et al. [28] considered rewards for the pairings that matched the crews' preferences in the objective function. ey solved the model via column generation and implemented it to a North American airline. e findings demonstrated that the suggested model could improve CRP solutions compared to previous studies. Mirjafari et al. [29] presented an integrated model for aircraft routing and CSP with the aim of minimizing crew, aircraft replacement, maintenance, and deadhead flight costs. ey suggest a new frame for maintenance operations based on flight hours. Because the Lagrangian relaxation method results have a lower gap to an optimum solution, this method is better than the particle swarm optimization (PSO) algorithm. A new integrated model for fleet assignment and crew scheduling problems was proposed by Rashidi Komijan et al. [1]. ey considered closed routes for crews and fleets simultaneously. To compare the effectiveness of the two methods, PSO and vibration damping optimization (VDO) algorithm solved ten large-scale examples. e findings showed that VDO leads to optimal results in a reasonable time. Saemi et al. [30] suggested an integrated CSP to minimize crew costs. e approach presented a better answer than computing CPP and CRP consecutively. It was also possible to add a day(s) off in a pairing to allow crew members to attend training courses and complete health checks. A multiobjective personalized model for the airline multiskilled was proposed by Baradaran and Hosseinian [31]. e multiskilled crew used in this model was assigned to flights and aircraft due to the crew's experience. ey used multiobjective differential evolution (MODE) and nondominated sorting genetic algorithm II (NSGA-II) to analyze the model. According to the results, the suggested model and algorithms could generate adequate schedules for CS problems, according to comparisons between algorithms. Ben Ahmed et al. [32] assigned the aircraft and crew to each flight simultaneously and included maintenance constraints in their model. ey solved the robust mixed-integer programming model using the decomposition approach and achieved good results. Table 1 briefly compares the proposed CRP model with the previous research studies.
To our best knowledge and Table 1, previous research studies mainly focused on cost minimization, fair workload distribution, and use of bidline and personalized approaches for crew rostering. is paper includes two main contributions. e first one is related to the issue of inconsistent (co)-pilots. is critical qualitative component should be evaluated to ensure passenger safety. e model prevents to use inconsistent crew in a single flight. e second contribution is that instead of suggesting preferred pairings, a scoring method is used. In contrast to previous studies, this strategy takes into account both the crew's willingness and unwillingness to assign to each pairing. In other words, one positive score is calculated for each favorable allocation and one negative score is considered for each undesirable allocation in the objective function. Total weighted preference is maximized by considering crew seniority.

Problem Statement
In this paper, a new mathematical model for the CRP is proposed. e objective function calculates the total weighted preference of the crew. We assume that each crew member declares his/her preferred and undesirable pairings.

Computational Intelligence and Neuroscience
By assigning a crew member to his/her preferred pairing, a positive score is calculated in the objective function. Similarly, if a crew member is assigned to his/her undesirable pairing, a negative score (penalty) is calculated. e objective function maximizes total weighted preference by considering crew seniority.
In each airline, some of the crew members may have bad relations with each other. Assigning inconsistent crew to a single flight may threaten flight safety. To cover this qualitative factor, the model should prevent to assign inconsistent crew to a flight. Also, several rules and regulations which are common in airlines are considered in the model. Some of these rules are as follows: (i) Time away from the base should not exceed a predefined value (ii) Maximum and minimum flying hours of crew are restricted (iii) An inexperienced co-pilot must be accompanied by an experienced pilot (iv) Some of the crew members should attend at training courses Due to the uncertainty of time away from the base (TAFB), the approach of Bertsimas and Sim [10] is used to formulate the robust counterpart. e advantage of using the approach of Bertsimas and Sim [10] is that the robust model will be linear. Also, the most well-known Indian airline, Air score cp , penalty cp ≥ 0.
e objective function maximizes total weighted preference of the crew. e way of calculating score and penalty is described in the constraints. Equations (2) and (3) ensure that each pairing is handled by one pilot and one co-pilot. Constraint (4) prevents assignment of inconsistent crew to a pairing. According to (5), the TAFB for each crew should not exceed the maximum allowed value. Relation (6) indicates that the total flying hours for a crew member should fall in the defined interval. Some pairings cannot be handled sequentially due to some reasons such as the destination of the first one is not the same as the origin of the next one. Constraint (7) prevents the assignment of a crew member to such pairings. Constraint (8) states that in each pairing, an inexperienced copilot must be accompanied with an experienced pilot. Relation (9) is about attendance of crew at the training course or meetings. Constraint (10) prevents concurrence of attendance at the training course and handling a pairing. Constraints (11) and (12) show the way of calculating score and penalty, which have been used in objective function. According to (11), one score is calculated if a crew member is assigned to his/her preferred pairing. Similarly, if he/she is assigned to an undesirable pairing, a penalty is considered. Constraint (13) ensures that a crew member can only be assigned to the pairings that origin from his/her home base. Finally, Constraints (14) and (15) indicate decision variable types. e seniority weight of pilots is in the range [0, 1]. Since the nature of this parameter is fuzzy, the Jimenez et al. [33] method has used to convert this uncertain coefficient in the objective function to definite. For more information, refer the study by Jimenez et al. [33]. For this matter, the triangular fuzzy distribution has chosen because of its simplicity and efficiency. After performing the calculations, the following equation replaces the indefinite objective function (1).
As TAFB is an uncertain parameter, the robust counterpart of the nominal model is formulated using the approach of Bertsimas and Sim [10]. TAFB includes total flying time, sit times, and layovers. TAFB p and TAFB p have assumed the nominal value and tolerance of time away from the base for pairing p. According to Bertsimas and Sim [10], protection level (Γ) is defined to make the model robustness. In other words, the solution's robustness is regulated by defining a protection level, which specifies a maximum deviation from the nominal TAFB. e main purpose of this parameter is to make the solutions generated workable for all uncertain scenarios. If at most Γ i of technical coefficients of the constraint i change in the defined interval, the model will definitely be robust. Also, a technical coefficient can be changed up to(Γ i − Γ i without violating the relation. Moreover, σ cp and θ cp are variables defined to formulate robust counterpart. To formulate robust counterpart, constraint (5) is replaced by the following constraints.

Solution Approach
According to NP-hard nature of the crew rostering problem, an appropriate optimization strategy is required [34]. Many rules and policies may deploy for crew rostering with Computational Intelligence and Neuroscience thousands of variables in a real-world problem. Furthermore, the presence of binary variables makes the solution approach more complicated.
Various papers in the literature use genetic algorithm to solve the airline crew scheduling problem. e GA is the second commonly applied solution technique and first metaheuristic approach for the CSP after column generation [35]. One of the benefits of the GA is the ability to manage numerous solution search areas [36]. As a consequence, the genetic algorithm is deployed as solution strategy for this study. In this section, problems are solved by a genetic algorithm in different sizes to evaluate the efficiency of the proposed approach. e GA outcomes are assessed and compared to the CPLEX solver in GAMS. Also, the solution method presented in this work and the components of the algorithm utilized are highlighted in Section 4.1.

Presentation of Solution.
e genetic algorithm uses the concept of chromosome to set parameters that offers a proposed problem solution. In practice, the chromosome is a string of solutions to solve the problem. e proposed chromosome can be designed as a string of discrete variables, binary numbers, and a continuous vector based on type of the problem. As a result, choosing an appropriate display form of the chromosome is an essential aspect of algorithm design. e designed chromosome for the proposed CRP is made up of the two-row vector. e number of columns in this vector is equal to the number of pairings |P|. e first and second rows specify the pilot number and the co-pilot number, respectively. Pilot Co-pilot Figure 1: Chromosome with 5 co-pilots and 6 pairings.

Signal-to-noise: Smaller is better
Maxlt npop pc pm

Generating the Initial Solution.
In this paper, the copilots are randomly assigned to each pairing. erefore, each compatible co-pilot gets selected to operate in pairing P. As an example, 6 pairings, 5 pilots, and 5 co-pilots are available (pilots are numbered from 1-5, and co-pilots are numbered from 6-10). An initial solution for the co-pilots assigned to the pairing is shown in Figure 1. To avoid assigning of a crew member to nonconsecutive pairings, Constraint (7) is applied to prevent the generation of infeasible solutions.
Pilot assignment process begins after all co-pilots get placed in their pairings. Each eligible pilot may get assigned to pairing P randomly, if the following conditions met: (1) e pilot should not be in conflict with the co-pilot (2) An inexperienced co-pilot must be accompanied by an experienced pilot in each pairing (3) Set of flights which may be assigned to a specific pilot in a sequence e pilot row is determined according to Figure 2 for the chromosome.

Generating a New Solution.
A suitable operator for generating a random neighborhood solution depends on to the nature of the problem. e CRP is an allocation problem; therefore, we considered the appropriate operators to generate new solutions. ese operators are defined in the following section.
(1) Mutation Operator. By using the swap operator, two pilots or two co-pilots are selected from two pairings and get swapped if the following conditions are met: (i) e co-pilot and the pilot are inconsistent (ii) If an inexperienced co-pilot is not paired with an experienced pilot (iii) Pairings that may not be done sequentially by one pilot or co-pilot For the mutation operator, a random number between 0 and 1 is generated. If the number is less than 0.5, a swap operation is performed on the co-pilots, and if it is greater than 0.5, it is performed on the pilots. Swapping the co-pilots of pairings P 5 and P 3 is presented in Figure 3.
(2) Crossover Operator. Figure 4 shows the single crossover operator for the CRP. Selecting the pilot/co-pilot row is similar to Section 4.1.3.1 for crossover operator as well.
e proposed chromosome in this paper complies with all constraints of the model except for the maximum/ minimum flight hours for crews and maximum time away from the base for crews. To tackle this problem, a penalty function strategy is applied. erefore, the average violation 1400 1200 100 -C 10 1400 1200 100 - Table 5: Crew assignment to pairings.

Pairings Crews Pilot
Co-pilot Computational Intelligence and Neuroscience 7 of these constraints is multiplied by a big number and then subtracted from the objective function.

Genetic Parameter Tuning.
In the genetic algorithm, maximum iterations, population size, crossover rate, and mutation rate parameters affect the performance of the algorithm. In this study, the Taguchi [37] is applied to get the best value for mentioned parameters. Each of the GA parameters is considered at three levels, as presented in Table 2.
Not all but only part of the factors may reach to their best level in Taguchi's [37] method. Each obtained result from Taguchi [37] method experiment is converting to a "signalto-noise" ratio. At this rate, the optimal value (average) is called the signal, and the undesirable value (standard deviation) is called the noise. According to the objective function which is to maximize the crew's satisfaction, the higher S/N means more desirability of results. erefore, the maximum point of the graph is selected for each parameter. Also, the corresponding level is considered as the optimal  e estimated values for GA parameters using the Taguchi [37] method are illustrated in Table 3 and Figure 5.

Computational Results for Test Problems.
e model is solved for a real case: Air India airline. Before solving the case study, 30 test problems in different scales are generated and solved using GAMS and GA to justify GA's efficiency in solving the proposed model. is operation is performed on a laptop with 16 GB RAM and Windows 10 64 bit. To illustrate the model results, a problem with 10 crew members is explained in detail. In this example, C 1 ′ , C 2 ′ , and C 3 ′ are assumed experienced pilots, while C 6 ′ and C 7 ′ are inexperienced co-pilots. Table 4 shows the data for this test problem. By solving the model, crew assignment is done. According to Table 5, one pilot and one co-pilot who are not inconsistent are assigned to each pairing. Crew's rosters can be easily concluded from Table 5. To illustrate, Figure 6 displays the first pilot's monthly roster. As the planning horizon is one    Table 6. In solving numerical examples by GAMS, the CPU time has been set to 10,800 seconds. According to Table 5, the performance of the GA to solve the proposed model is justified due to two reasons: (1) the maximum optimality gap between the genetic algorithm and exact solution is only 1.482%, and (2) in medium and large-sized problems, the CPU times in the GA are considerably less than GAMS.
A comparison between the CPU times in GA and GAMS is shown in Figure 7. As maximum solution time has been set to 10,800 seconds in GAMS, computations in problems 14-30 are terminated after three hours.
A comparison of the results between the results of GAMS and genetic algorithm is shown in Figure 8. e results of the objective function are shown in Figure 8. e green and orange points represent the objective function values resulted from GA and GAMS, respectively. It is worth noting that GAMS solution time has been limited to three hours. As a result, GAMS did not reach the optimum solution in problems 14 to 30.

Case Study.
Air India is a well-known airline in India, headquartered in New Delhi. As one of the biggest international airlines, Air India has 11% of the domestic market share [38]. Air India's main hub is Indira Gandhi International Airport in New Delhi, with a secondary hub at Chhatrapati Shivaji Maharaj International Airport in Mumbai. According to reports in 2019, Air India fleets include Airbus and Boeing that serve 170 destinations all over the world [32].
To solve the proposed model for the case study, a genetic algorithm was applied and the Taguchi [37] approach was used to adjust its parameters.
is case included 6190 pairings and 1340 (co)-pilots. e model was solved in 761.188 seconds, and the optimum value of the objective function was 5034.002. e convergence of the GA results is shown in Figure 9. Converging after 400 iterations indicates good solution quality.
As discussed in Section 3, the objective function includes two terms: score and penalty for assigning crew to desirable/ undesirable pairings. According to the Air India results, crew members were assigned to their preferred pairings in 72% of assignments. Table 7 shows the monthly rosters of 20 experienced pilots, and 20 co-pilots as a small sample. Moreover, the attendance days of 90 co-pilots in the required training courses are presented in Table 8.
To compare the performance of the proposed model with the previous research studies in which undesirable assignment penalty had not been considered, penalty term was removed from the objective function and constraint (12) was also disabled. By removing penalty concept, 2065 undesirable pairings are assigned to crews. is means that each crew will have an average of at least two undesirable pairings per month. However, our proposed model has reduced the number of undesirable rosters to 739, which is about 61.59% less than the previous models.

Robust Optimization.
In this section, Bertsimas and Sim [10] approach is used to formulate the robust counterpart of the proposed model. e approach of Bertsimas and Sim leads to a linear robust model, while it is independent of knowing the distribution of data and control of the level of protection through parameter Γ [39].
TAFB is subject to several factors such as flight time, flight duration, and connection time, which makes it an uncertain parameter. By considering different values for Γ (0 ≤ Γ ≤ 1), the robust model was solved for numerical

Computational Intelligence and Neuroscience
′ × example 3, and the results of GAMS and GA were compared. After justifying the performance of the GA, the robust model was applied to the case study and the results were analyzed ( Table 9).
As represented in Table 9, the optimality gap of GA for numerical example 3 is negligible. erefore, the performance of GA to solve the robust model is acceptable. According to the robust model applied to the case study, the model acts more conservatively by increasing of the Γ. In the case of Γ � 1, it guarantees that feasible solution to be at its maximum value and the objective function to be at its minimum value. Figure 10 shows the total weighted desirability obtained from the genetic algorithm for 400 iterations.

Conclusions
Crew scheduling problem is often divided into CPP and CRP and solved consecutively to reduce complexity. In this article, a novel MIP model was introduced for CRP. e difference between the proposed model and previous research is considering a scoring mechanism instead of bidding for preferred pairings, crew seniorities, and inconsistent crew. In addition, the model considers common rules and regulations as well as crew attendance at the required training courses. With this advantage, obviously, better solutions are obtained. Moreover, the robust counterpart of the nominal CRP model was presented to the uncertainty of mission time of the crews (TAFB). In addition, the suggested model's efficiency was assessed using real-world data from Air India Airline. Because the model was NP-hard, the CPLEX solver in GAMS was used for 13 first test problems, and the GA was utilized for 30 test problems with different scales and the case study. e results show that GA has a 1.482% average gap in optimum solution compared to exact solution. e proposed CRP is computed utilizing Air India data, which is one of the largest international carriers. Air India provides services to 170 destinations in 31 countries. e proposed model was tested for 6190 pairings and 1340 crews. For this real-world data, the objective function value was 5034.002 with 761.188 seconds of computational time. is research has two limitations that can be improved in future research studies. e first one is that (co)-pilots' medical checks have not been considered. e second one is that reserve crew have not been included in the model. Moreover, the following suggestions can be considered for future research [40][41][42]: (i) Considering a fair distribution of working time among crew members (ii) Providing an integrated model for CRP with other steps of the airline scheduling problem, and considering disruption for the crew rostering (iii) Developing the cooperative game between airlines to exchange crews among them

P:
Set of all pairings p: Index of pairing (p � 1, 2, . . .) C: Set of all crew members (pilots and co-pilots) c, c ′ : Index of pilots and co-pilots (c, c ′ �1, 2, . . .) C 1 : Set of all pilots c 1 : Index of pilots (c 1 � 1, 2, . . ., C 1 ) C 2 : Set of all co-pilots c 2 : Index of pilots (c 2 � 1, 2, . . ., C 2 ) C 1 ′ : Set of all experienced pilots c 1 ′ : Index of experienced pilots (c 1 ′ �1, 2, . . ., C 1 ′ ) C 2 ′ : Set of all inexperienced co-pilots c 2 ′ : Index of inexperienced co-pilots (c 2 ′ � 1, 2, . . ., C 2 ′ ) C inc : Set of inconsistent crew C inst : Set of crew who must participate at the training course P d : Set of all pairings that include day d NA: Set of all pairings that cannot be performed sequentially by a (co)-pilot D: Set of all days in which the training course is held d: Index of training course (d �1, 2, . . .) Binary parameter that equals 1 if pairing p is desirable for (co)-pilot c, and 0 otherwise Base cp :

Parameters
Binary parameter that equals 1 if pairing p begins and ends at the home base of (co)-pilot c, and 0 otherwise w c : Seniority weight of (co)-pilot c Decision variables X cp : Binary variable that equals 1 if pairing p is assigned to (co)-pilot c, and 0 otherwise Y cd : Binary variable that equals 1 if (co)-pilot c attends at training course held on day d, and 0 otherwise score cp : Score of assigning (co)-pilot c to his/her preferred pairing p penalty cp : Penalty for assigning (co)-pilot c to the undesirable pairing p. Computational Intelligence and Neuroscience 13 Data Availability e data used to support this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.