Patrick A Pintus
email: patrick.pintus@univ-amu.fr.

Risk-Taking, Global Diversification and Growth: Comment *

Keywords: International Financial Integration, Endogenous Growth, Stochastic Stability Journal of Economic Literature Classification Numbers: F34, F43, O40

 showed that growth and welfare gains from international risk-sharing arise in a continuous-time stochastic AK model. More precisely, he proved that a portfolio shift from safe and low-return capital to riskier and high-return capital triggers an unambiguous increase in growth. In this note I stress necessary and sufficient conditions ensuring stochastic stability of the exponential balanced-growth path, an issue that has not been addressed by Obstfeld. Not surprisingly, stability requires the average of the wealth growth rate to be positive, which makes clear how mean growth should be defined. Differently, though, Obstfeld defines mean growth as the growth rate of average wealth, which is larger than the mean growth rate of wealth under the maintained assumption that wealth is log-normally distributed, because the latter growth concept is risk-adjusted. The two notions of mean growth have very different comparative statics properties both for economies that hold some risk-free capital and for economies that fully specialize in risky capital. Different from Obstfeld's results, international financial integration increases the stability-related mean growth rate for both complete and incomplete specialization, if risk aversion takes on moderate values and provided that the intertemporal substitution elasticity is smaller than one. Although the welfare computations presented by Obstfeld are preserved, because they ultimately depend on parameter values, this note shows that stochastic stability sheds new light on the mechanisms that trigger growth changes under financial integration and underlines the intuition behind them.

Introduction

Two decades ago was published a seminal article by [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF], which shows how large welfare gains from international risk-sharing arise in a continuous-time stochastic AK model. The intuition for why such large welfare benefits materialize is developed first, in section I of that paper, in the closed-economy version of the model where individuals may invest in risk-free capital and risky capital. The central result is that if the economy holds some risk-free capital, a fall in exogenous risk unambiguously leads to an increase in the share of wealth invested in risky capital. In other words, a portfolio shift toward riskier capital triggers a return effect that dominates the possibly adverse effect of a higher propensity to consume out of total wealth. As a consequence, growth and welfare go up. In contrast, goes the intuition, complete specialization in risky capital results in an ambiguous effect of exogenous risk on growth, that is entirely governed by whether the intertemporal substitution elasticity is larger or smaller than one.

In this note I show that such results are based on a definition of mean growth that seems natural absent any concern about stability of the balanced-growth path but that turns out to be misleading if the issue of stability is addressed, as it should be. I make clear why well-known conditions ensuring stochastic stability of the balanced-growth path (exposed in [START_REF] Mao | Stochastic differential equations and applications[END_REF]used in Boucekkine and[START_REF] Boucekkine | Stochatic stability of endogenous growth: the AK case[END_REF] relate to an alternative definition of mean growth. To understand why a confusion may arise, it is important to stress that [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] considers an AK model where risk is introduced via diffusion equations that involve geometric Brownian motions. In such a setting, the dynamics of the economy that result from optimal portfolio and savings decisions are described by a linear, stochastic differential equation with fixed coefficients and one might wrongly think that two correct concepts of growth are available to freely choose from. The first defines mean growth as the growth rate of average wealth, which is the one used by [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] and many others (see e.g. [START_REF] Jones | Stochastic models of growth[END_REF]Manuelli, 2003, Steger, 2005). The alternative definition, which emerges naturally in the quest for stability conditions that is developed below, is to define mean growth as the average growth rate of wealth, as opposed to the growth rate of average wealth. Because by assumption wealth is lognormally distributed, it follows by Jensen's inequality that the average growth rate of wealth -second notion -is lower than the growth rate of average wealth -first notion.

Not surprisingly, stochastic stability holds if and only if the average growth rate is positive, a condition that is stronger than the requirement that the growth rate of average wealth be positive.

The contribution of this note is to show that very different comparative statics results obtain when one uses the second definition of mean growth, as one should in view of stability conditions. More precisely, mean growth happens to be enhanced by financial integration under conditions that would possibly lead to the opposite conclusion if one were to use the definition of mean growth advocated in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF]. This property is most striking in a specialized economy, where for example a fall in exogenous risk results in larger growth even if the intertemporal substitution elasticity is smaller than one, despite the fact that a portfolio shift does not happen. Before moving on to the technical analysis, I should stress that the welfare computations provided by [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] are completely unaltered by stability considerations, because they of course ultimately depend on parameter values and not on how one defines mean growth. However, stability unequivocally imposes a definition of mean growth that leads to very different comparative statics and, therefore, sheds new light on the mechanisms that trigger growth changes. For instance, it turns out that both risk aversion and the intertemporal substitution elasticity matter to determine whether or not financial integration boosts growth. Section 3, in particular, shows how conflicting return and variance effects help explain what happens in intuitive terms. A summary of the results unveiled in this note is that although growth gains from international financial integration turn out to be smaller than previously thought, they are likely to be positive both for specialized economies and for economies that hold some risk-free capital, under reasonable parameter values.

The paper is organized as follows. Section 2 presents conditions for stochastic stability of the exponential balanced-growth path and the associated definition of mean growth.

Section 3 then performs comparative statics analysis and compares results with those obtained by [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF].

2 Stochastic Stability and the Definition of Mean Growth [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] considers an AK version of the optimal portfolio model developed in [START_REF] Merton | Lifetime portfolio selection under uncertainty: the continuous-time case[END_REF] and [START_REF] Samuelson | Lifetime portfolio selection by dynamic stochastic programming[END_REF]. The following equation describes optimal wealth accumulation and is identical to equation [14] derived in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF]

1 : dW = [ωα + (1 -ω)i -μ]W dt + ωσW dz (1)
where α > 0 and i > 0 are the mean returns of risky capital and risk-free bonds respectively, μ > 0 is the average propensity to consume out of total wealth, dz is a Wiener process and σ 2 ≥ 0 is the exogenous variance of the return on risky capital . In equation ( 1), 1 ≥ ω ≥ 0 denotes the share of wealth invested in risky capital and its expression is given in equation [11], that is:

ω ≡ α -i Rσ 2 > 0 ( 2 )
Two cases occur, depending on whether ω is smaller than or equal to one. We will refer to the first case as incomplete specialization -when the economy holds some risk-free bond -and to the second case as complete specialization -when the economy has all its wealth invested in risky capital. It is important to notice a major difference between the two configurations: when ω < 1, [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] shows that i equals r, the mean return on risk-free capital such that r < α, so that a fall in exogenous risk σ 2 always results in a portfolio shift away from risk-free capital, that is, ω goes up. This first case occurs when Rσ 2 > αr, that is if (utility adjusted) risk is large enough to prevent complete specialization. If, however, Rσ 2 < αr, specialization is complete because risk is small enough to ensure ω = 1. In that case a fall in exogenous risk triggers a rise in risk-free return i = α -Rσ 2 that compensates for the fall in σ 2 so that the economy keeps all its wealth in risky capital and enjoys lower risk.

Conditions for stochastic stability of linear differential equations with fixed coefficients are stated in chapter 4 of [START_REF] Mao | Stochastic differential equations and applications[END_REF] and applied in [START_REF] Boucekkine | Stochatic stability of endogenous growth: the AK case[END_REF] to a simpler AK model. Stated in terms of a generic diffusion equation dX = φXdt + θXdz, the unique solution of which is X(t) = X(0) exp{(φ -0.5θ 2 )t + θz(t)}, the main result is that the variable X tends exponentially to infinity with probability one, as time goes to infinity, if and only if φ -0.5θ 2 > 0, whereas X tends to zero with probability one if and only if φ -0.5θ 2 < 0. A straightforward application of such conditions leads to the following lemma:

Lemma 2.1 (Stochastic Stability of the Balanced-Growth Path)

Wealth tends exponentially to infinity, along a balanced growth path, with probability one when time tends to infinity if and only if ωα

+ (1 -ω)i -μ > 0.5ω 2 σ 2 .
Obstfeld (1994) defines mean growth -g in his notation -as the growth rate of average wealth, which is given in view of equation ( 1) by ωα + (1ω)iμ or, equiva-lently after plugging the expression of the share invested in risky capital, by g(

E [W ]) = ε(i-δ)+(1+ε)(α-i) 2 /(2Rσ 2 )
, where ε > 0 is the elasticity of intertemporal substitution in consumption and R > 0 is relative risk aversion (see equation [16] in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF].

Lemma 2.1 shows that g(E[W ]) > 0 is necessary but not sufficient for stochastic stability of exponential growth at a positive rate. In other words, assuming g( E [W ]) > 0 would result in convergence to zero wealth with probability one provided that g(

E [W ]) < 0.5ω 2 σ 2 .
This is an example of the well-known fact that noise, if big enough, can significantly alter and sometimes overturn convergence: although in that case the deterministic model experiences exponential growth at a positive rate, the stochastic version converges to the trivial solution W = 0 with probability one (see [START_REF] Mao | Stochastic differential equations and applications[END_REF][START_REF] Boucekkine | Stochatic stability of endogenous growth: the AK case[END_REF][START_REF] Boucekkine | Stochatic stability of endogenous growth: the AK case[END_REF], for examples of that sort).

Lemma 2.1 therefore suggests that mean growth should be defined as the average of the wealth growth rate2 , that is, E[g(W )] ≡ ωα + (1ω)iμ -0.5ω 2 σ 2 which can be simplified, using the expression of ω in (2), to:

E[g(W )] = ε(i -δ) + (α -i) 2 2Rσ 2 1 + ε - 1 R (3)
where δ > 0 is the subjective rate of time preference. A few comments are in order.

Because wealth is assumed to be log-normally distributed, the property that

E[g(W )] < g(E[W ]
) follows, of course, from Jensen's inequality: the expected value of the log of wealth is smaller than the log of expected wealth and a similar inequality applies to their derivatives with respect to time. More importantly, one goes from the first definition of mean growth, used by [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF], to the second, more appropriate, one by subtracting half of the (endogenous) variance of wealth, that is, (αi) 2 /(2R 2 σ 2 ), hence the additional term -1/R in equation (3). Therefore, comparative statics results are expected to be very different, as I show next.

Comparative Statics of Mean Growth

Straightforward computations lead to the following main result of this note.

Proposition 3.1 (Comparative Statics of Mean Growth)

The dynamics of wealth accumulation defined in equation ( 1) has two regimes: Not surprisingly, comparing Proposition 3.1 and results in Obstfeld (p. 1315Obstfeld (p. , 1994) ) shows important differences. As shown in case (i) of Proposition 3.1, incomplete specialization results in larger growth when exogenous risk falls down only for large enough values of either risk aversion or of the intertemporal substitution elasticity. In contrast, [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] claims that a portfolio shift unambiguously improves growth, independent of R and ε. When the correct definition of mean growth is used, this is no longer true. In addition, case (ii) of Proposition 3.1 shows that the results obtained by [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] for complete specialization can be overturned under reasonable assumptions on parameters. A striking example is the case of unitary intertemporal substitution elasticity, that is, ε = 1. Whereas this case implies that the growth rate of average wealth is independent of exogenous risk in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] (see his equation [17]), case (ii) in Proposition 3.1 shows that the average growth rate is in fact a decreasing function of exogenous risk for all values of risk aversion. This property suggests that international financial integration is likely to boost growth in economies that invest all their wealth in risky capital.

(i) if Rσ 2 > α -r (incomplete
More generally, conditions ensuring that a fall in exogenous risk boosts growth for both complete and incomplete specialization become clearer under the assumption that the intertemporal substitution elasticity is smaller than one, which seems to accord better with empirical measures. Remember that in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF], in this case growth unambiguously goes up under incomplete specialization whereas growth slows down in specialized economies, following financial integration. In contrast, Proposition 3.1 shows that using the correct definition of mean growth delivers a more contrasted picture: when ε < 1, a fall in exogenous risk leads to an increase in mean growth provided that relative risk aversion takes on moderate values, that is, if and only if 1/(1ε) > R > 1/(1 + ε). For example, the latter inequalities are met when R = 1. The bottom line is that because it leads to smaller exogenous risk, financial integration is expected to improve mean growth for both complete and incomplete specialization under reasonable parameter values.

So as to clarify the intuition behind the striking differences with results reported in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF], I now focus on the case such that ε = 1, which leads to the well-known result that the average propensity to consume out of total wealth is then given by the impatience rate, that is, μ = δ. This assumption neutralizes the effect of exogenous risk on the consumption-wealth ratio, which has been described in earlier papers and in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] in particular. I now explain how a fall in exogenous risk affects mean growth. Again, two cases arise depending on the level of exogenous risk.

(i) if Rσ 2 > αr, specialization is incomplete because exogenous risk is so large that the economy holds some risk-free capital (that is, ω < 1). It follows that the risk-free interest rate i = r and that the expression for mean growth simplifies to:

E[g(W )] = r -δ + (α -i) 2 Rσ 2 mean return effect -0.5 (α -i) 2 R 2 σ 2 variance effect (4)
The expression for mean growth in equation ( 4) reveals that two conflicting effects are at work. The return effect is such that a fall in exogenous risk σ 2 boosts welfare growth because, as shown by [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF], the portfolio shift away from risk-free capital increases growth under the assumption that risky capital has a larger mean return than risk-free capital. However, although ignored by [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF], a variance effect also materializes, essentially because a larger share in the risky asset implies that the endogenous variance of wealth goes up when exogenous risk goes down. Stochastic stability of the balanced-growth path requires the variance effect to be not too large but such a condition does not exclude that mean growth be a decreasing function of exogenous risk, then overturning Obstfeld's result, if risk aversion is less than one half.

(ii) if αr > Rσ 2 , specialization is complete (ω = 1). It follows that the risk-free interest rate adjusts to ensure i = α -Rσ 2 > r and that the expression for mean growth simplifies to:

E [g(W )] = α -δ - 0.5σ 2 variance effect (5)
Equation ( 5) makes clear what happens when specialization is complete. In contrast to case (i), there is no return effect because the economy already benefits from full specialization so that a fall in σ 2 has no effect on the mean return -there is no portfolio shift. However, a variance effect still occurs but it now has an opposite effect on mean growth compared to case (i). This is because the endogenous variance of wealth now goes down when exogenous risk goes down, as the risk-free return goes up to ensure that specialization remains complete in the face of a fall in risk. Quite interestingly, an analysis based on the alternative but misleading notion of mean growth, as in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF], predicts that growth is not affected by such a fall in risk.

Relaxing the assumption that ε = 1 delivers similar intuitions. In case (i) the return effect dominates the variance effect so that a fall in exogenous risk fosters growth if and only if risk aversion is large enough. In case (ii) there is no return effect and the variance effect, now working in opposite direction, implies that growth improves after a fall in exogenous risk only if risk aversion is not too large when the intertemporal substitution elasticity is smaller than unity. In other words, our results about specialized economies accord with the well-documented trade-off between growth and volatility under reasonable assumptions about attitudes toward risk, for example if relative risk aversion equals one. In contrast, incomplete specialization leads to a positive relationship between the mean growth and variance of wealth under unitary risk aversion. Overall these results

suggest that taking into account the variance effect on mean growth, which has been ignored by [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF], yields the prediction that the effects of financial integration on economies that specialize in risky capital do not qualitatively differ from those on economies that hold some risk-free capital if reasonable parameter values are assigned to risk aversion and intertemporal substitution.

To make the comparison even more transparent, I now reproduce and extend in Table 1 a numerical example given in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF]. More precisely, Table 1 starts with the Example 1 that is presented in pages 1318-1319 of [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] and that assumes R = 4 and ε = 0.5 in particular. Table 1 compares the magnitudes of both definitions of mean growth under this parameterization and also, for robustness purpose, when R = 1 while all other parameter values are unchanged.

( E [W ]) in Right Panel R = 4 R = 1 Autarky 1.41% 1.25% Integration 1.75% 1.38% R = 4 R = 1 Autarky 1.69% 1.75%
Integration 2.00% 1.63%

In line with the analytical characterization outlined above, comparing both panels in Table 1 confirms that the mean growth rate of wealth is lower than the growth rate of mean wealth. More interestingly, comparing the rightmost columns of both panels reveals that, when R = 1, the conclusion regarding growth that is obtained by [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] is overturned when the appropriate concept of mean growth is adopted. In fact, while the right panel predicts that growth falls (by about 12 basis points) after integration in the case of full specialization, it turns out that growth actually goes up (by about 13 basis points) as depicted in the left panel that uses the appropriate definition of mean growth. Let me stress that although welfare computations reported in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] are not altered at all by stability considerations, the examples in Table 1 further confirm that different comparative statics properties obtain when the stability-related concept of mean growth is used, as it should be. Aside from theoretical concerns, this is also relevant for empirical research, which typically aims at measuring the growth gains from international financial integration.

Conclusion

This note shows, by way of analytical results and numerical examples, that the comparative statics results derived in [START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] are misleading because they are based on an inappropriate notion of mean growth. Conditions ensuring that the exponential balanced-growth path is stable, in the stochastic sense, reveal that mean growth should be defined as the average growth rate of wealth, as opposed to the growth rate of average

wealth. With such a definition in hand, I show that international financial integration leads to very different comparative statics results and that it is much more likely to boost growth, both for fully specialized economies that invest all their wealth in risky capital and for economies that hold some risk-free capital. This note shows that although the growth gains from financial integration are in fact smaller than previously thought, they are more likely to be positive for empirically reasonable values of relative risk aversion and the intertemporal substitution elasticity, independently of whether or not the economy invests in risk-free capital.

  specialization): the average growth rate E [g(W )] is a decreasing function of exogenous risk σ 2 if and only if R(1 + ε) > 1, that is, for large values of either risk aversion or of the intertemporal substitution elasticity.(ii) if αr > Rσ 2 (complete specialization): the average growth rate E[g(W )] is a decreasing function of exogenous risk σ 2 if and only if R(1ε) < 1, that is, for small values of risk aversion and large values for the intertemporal substitution elasticity.

Table 1 .

 1 Numerical Values of E [g(W )] in Left Panel and g

For clarity, I use square brackets to label equations that appear in[START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] and round brackets for equations in this paper.

In fact,[START_REF] Obstfeld | Risk-taking, global diversification and growth[END_REF] uses later in his paper this notion for measurement purpose, e.g. in page 1321, although not for the comparative statics analysis developed at the outset.


Nathan Krislock
email: krislock@math.niu.edu

Jérôme Malick
email: jerome.malick@inria.fr

Frédéric Roupin
email: frederic.roupin@lipn.univ-paris13.fr

Computational results of a semidefinite branch-and-bound algorithm for k-cluster

Keywords: combinatorial optimization, semidefinite programming, triangle inequalities, k-cluster problem, k-densest subgraph problem

This computational paper presents a method to solve k-cluster problems exactly by intersecting semidefinite and polyhedral relaxations. Our algorithm uses a generic branch-and-bound method featuring an improved semidefinite bounding procedure. Extensive numerical experiments show that this algorithm outperforms the best known methods both in time and ability to solve large instances. For the first time, numerical results are reported for k-cluster problems on unstructured graphs with 160 vertices.

Introduction

The k-cluster problem

Given a graph G = (V, E), the k-cluster problem consists of determining a subset S ⊆ V of k vertices such that the sum of the weights of the edges between vertices in S is maximized. This is a classical problem of combinatorial optimization, also known under the names "heaviest k-subgraph problem", "k-dispersion problem", and "k-densest subgraph problem" (when all the weights are equal to one). The problem can be seen as a generalization of the max-clique problem, and also as a particular case of quadratic knapsack problem where all the costs are equal.

Letting n = |V | denote the number of vertices, and w ij denote the edge weight for ij ∈ E and w ij = 0 for ij / ∈ E, the problem can be modeled as the 0-1 quadratic optimization problem (KC) maximize

1 2 z T W z subject to n i=1 z i = k z ∈ {0, 1} n , ( 1
)
where W := (w ij ) ij is the weighted adjacency matrix of the graph G.

The k-cluster problem (1) is a fundamental graph optimization problem, and arises in many applications such as telecommunication, warehouse location, military defence, social networks, and molecular interaction networks; see more details and references in the introductions of, e.g., [Pis06, MR12, BCV + 12].

It is well-known that the k-cluster problem is a hard combinatorial optimization problem: it is NP-hard (even for special graphs, see,e.g., [START_REF] Feige | Approximation algorithms for maximization problems arising in graph partitioning[END_REF]), it does not admit a polynomial time approximation scheme [START_REF] Corneil | Clustering and domination in perfect graphs[END_REF][START_REF] Khot | Ruling out PTAS for graph min-bisection, dense ksubgraph, and bipartite clique[END_REF], and there is even a huge gap between the best known approximation algorithm and the known inapproximability results (see, e.g., [BCC + 10, BCV + 12]).

In practice, k-cluster problems are also difficult to solve to optimality. Even if there are many theoretical articles on k-cluster, there are only a few of them on exact resolution. Among the only works attacking this problem are the pioneering [START_REF] Erkut | The discrete p-dispersion problem[END_REF], LP-based branch-and-bound method of [Pis06, Sec 2.4], and the convex quadratic relaxation method of [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] using semidefinite programming and CPLEX. Note that, before 2006 and [START_REF] Pisinger | Upper bounds and exact algorithms for p-dispersion problems[END_REF], no non-trivial k-cluster problem of size n = 100 was able to be solved. As of 2014, the state-of-the-art method for solving problem (1) to optimality is the semidefinite-based branch-and-bound algorithm of [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF] which is able to solve instances of size n = 120, and which compared positively to the previous methods on smaller problems.

Contribution and outline of this article

Our recent work [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] presents an improved semidefinite bounding procedure for Max-Cut, another classical combinatorial optimization problem. Max-Cut problems can be written as maximizing a quadratic function over the vertices of a hypercube, that is, as (1) but without the equality constraint n i=1 z i = k.

In this current paper, we build upon both [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF] and [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] by adapting and extending for the k-cluster problem the bounding procedure of [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF]. As we will see, extending techniques to k-cluster that have proven effective for Max-Cut brings complications due to the presence of the additional linear constraint. However, the extensive numerical experiments of this paper show that the resulting algorithm greatly outperforms the methods of [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF][START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF], which are the previous best existing methods to solve the k-cluster problem to optimality. Our algorithm is also able to solve unstructured k-cluster problems of sizes n = 140 and n = 160, for which no numerical results have been reported in the literature. The main contribution of this paper is thus to advance our ability to solve k-cluster problems to n = 160 from the previous limit of n ≤ 120.

The outline of this paper is as follows. In Section 2 we describe our improved semidefinite bounding procedure for the k-cluster problem (1). In Section 3 we describe our branch-and-bound implementation using our improved semidefinite bounding procedure for solving k-cluster problems to optimality. In Section 4 we present our numerical results. Finally, we give concluding remarks in Section 5.

Improved semidefinite bounding procedure for k-cluster

In this section we describe the improved bounding procedure that is based on semidefinite programming (SDP) bounds of k-cluster. We start by briefly recalling the standard strengthened SDP relaxation of k-cluster that we will approximate with our bounds. We use the following standard notation: the inner product of two matrices X and Y is X, Y := trace(X T Y ), and X 0 means that X is symmetric positive semidefinite.

Strengthened semidefinite relaxation

Semidefinite relaxations of the k-cluster problem (1) have already been considered in many papers for different purposes, such as [START_REF] Roupin | From linear to semidefinite programming: an algorithm to obtain semidefinite relaxations for bivalent quadratic problems[END_REF] for a computational study of the bounds, [BCV + 12] for (in)approximation results, and [START_REF] Pisinger | Upper bounds and exact algorithms for p-dispersion problems[END_REF] and [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF] in the context of exact resolution.

The algorithm presented in this paper uses a strengthened semidefinite relaxation in {-1, 1} variables as in [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF]. The derivation of this SDP relaxation uses standard techniques (see, e.g., [START_REF] Poljak | A recipe for semidefinite relaxation for (0,1)-quadratic programming[END_REF][START_REF] Saigal | Handbook of Semidefinite Programming[END_REF]), namely reinforcement by redundant constraints, homogenization, change of variables, and lifting to the space of matrices -see the reformulations 1-3 in [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF] Section 1] leading to the SDP relaxation [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF]Equation (9)]. For the sake of brevity, we do not repeat here the derivation of the relaxation, and we describe only the final SDP problem, refering to [MR12, Section 2] for more modeling information.

We consider the bound for the k-cluster problem (1) given by the following SDP problem:

(SDP I ) maximize Q, X subject to Q j , X = 4k -2n, j ∈ {0, . . . , n}, diag(X) = e, X 0, A I (X) ≥ -e, (2)
where X lies in S n+1 , e is the vector of all ones, and

Q := 1 4 e T W e e T W W e W , Q j := 0 e T + (n -2k)e T j
e + (n -2k)e j e j e T + ee T j , for j ∈ {0, . . . , n}, with e j ∈ R n being the j-th column of the n × n identity matrix, for j ∈ {1, . . . , n}, and e 0 := 0 ∈ R n . To lighten notation, we will gather all the equality constraints (including the diagonal ones) together as B(X) = b, where b ∈ R 2n+2 and B : S n+1 → R 2n+2 is a linear operator. Let us explain briefly the role of the constraints; again, for more details about this formulation, we refer to [MR12, Section 1]. The j = 0 constraint in problem (2) comes from the constraint n i=1 z i = k in the original k-cluster problem (1). In addition, we further strengthen the SDP bound by adding reinforcing equality constraints and valid triangle inequality constraints: i) Reinforcing equality constraints. The reinforcing equality constraints,

Q j , X = 4k -2n, j ∈ {1, . . . , n} , (3)
come from including the redundant product constraints n i=1 z i z j = kz j , for j ∈ {1, . . . , n}, to the original problem before forming the SDP (Lagrangian) relaxation. It is known that the reinforcing equality constraints (3) provide the best possible SDP bound when considering the inclusion of valid equality constraints (see, e.g., [START_REF] Faye | Partial Lagrangian for general quadratic programming. 4'OR[END_REF]). More precisely, when adding any set of redundant quadratic constraints {z T C j z + b T j z + a j = 0 : j ∈ J} before forming the SDP (Lagrangian) relaxation, the resulting bound is greater or equal to the partial Lagrangian relaxation of k-cluster,

(DP) min µ max z s.t. e T z=k 1 2 z T W z + µ i (z 2 i -z i ) ,
where only the binary constraints z ∈ {0, 1} n (written equivalently as z 2 iz i = 0, for i = 1, . . . , n) are relaxed, but the equality constraint n i=1 z i = k (written equivalently as e T z = k) is not relaxed. It is shown in [START_REF] Faye | Partial Lagrangian for general quadratic programming. 4'OR[END_REF] that the set of product constraints achieves the partial Lagrangian relaxation bound, and is therefore optimal when considering the inclusion of valid equality constraints. ii) Triangle inequalities. The triangle inequalities are defined by

X ij + X ik + X jk ≥ -1, X ij -X ik -X jk ≥ -1, -X ij + X ik -X jk ≥ -1, -X ij -X ik + X jk ≥ -1,
for 1 ≤ i < j < k ≤ n + 1, and correspond to the fact that for any x ∈ {-1, 1} n+1 , it is not possible to have exactly one of three products {x i x j , x i x k , x j x k } equal to -1, nor is it possible to have all three of the products equal to -1. There are a large number, 4 n+1 3 , of triangle inequalities in total. Therefore, we will iteratively add a subset of the most violated inequalities. For a subset of triangle inequalities I, we let A I : S n → R |I| be the corresponding linear function describing the inequalities in this subset. For every set of triangle inequalities I, the SDP relaxation in problem (2) gives us an upper bound on the value of the maximum weight k-cluster:

(KC) ≤ (SDP I ). ( 4
)
As was shown in [START_REF] Roupin | From linear to semidefinite programming: an algorithm to obtain semidefinite relaxations for bivalent quadratic problems[END_REF], the (SDP I ) bound with all triangle inequalities (i.e., |I| = 4 n+1

3 ) is tight in that it achieves the optimal value of the k-cluster problem on some instances, but it is also expensive to compute. Instead of using the (SDP I ) bound directly, our approach here is based upon the semidefinite bounds of [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds for 0-1 quadratic problems leading to quasi-newton methods[END_REF] that we present in the next section.

Before moving on to present our semidefinite bounds, we point out an easy but important result on the SDP relaxation (2): unlike the SDP relaxation of Max-Cut used in [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF], problem (2) is not strictly feasible.

Lemma 1. The semidefinite relaxation (2) of the k-cluster problem is not strictly feasible (i.e. a feasible matrix cannot be positive definite).

Proof 1. Let X be feasible for problem (2). Let v := 4k -2n -2e ∈ R n+1 . By partitioning X consistently with the vector v, we have

Xv = X 11 X T 21 X 21 X 22 4k -2n -2e = (4k -2n)X 11 -2X T 21 e (4k -2n)X 21 -2X 22 e .
Since diag(X) = e, we have X 11 = 1. The j = 0 constraint can be written as

0 e T e 0 , X 11 X T 21 X 21 X 22 = 4k -2n.
Simplifying we get 2X T 21 e = 4k -2n. Thus, the first component of Xv is zero:

(4k -2n) X 11 -2X T 21 e = (4k -2n) -(4k -2n) = 0.
Next we let j ∈ {1, . . . , n} and consider constraint j:

0 e T + (n -2k)e T j e + (n -2k)e j e j e T + ee T j , X 11 X T 21 X 21 X 22 = 4k -2n.
Simplifying we get 2X T 21 e -(4k -2n)e T j X 21 + 2e T j X 22 e = 4k -2n. Since 2X T 21 e = 4k -2n, we have (4k -2n)e T j X 21 -2e T j X 22 e = 0. Since this holds for all j ∈ {1, . . . , n}, the last n components of Xv are also zero:

(4k -2n)X 21 -2X 22 e = 0.
Thus we can conclude that Xv = 0. Since v = 0, we have that X is not positive definite. Therefore, problem (2) is not strictly feasible.

In conic programming, strict feasibility is often used to guarantee no duality gap; see, e.g., [BTN01, Th.2.4.1]. In spite of the previous lemma, there is still no gap between (2) and its dual

minimize b T y + e T z subject to Q -B * (y) + A * I (z) 0, z ≥ 0, (5)
where B * and A * I are the adjoint operators of B and A I respectively. Here it is the strict feasibility of the dual problem that implies zero duality gap. Since we have for y ∈ R 2(n+1)

B * (y) = n+1 j=1 y j Q j-1 + Diag(y n+1 , . . . , y 2n+2 ) (∈ S n+1 ),
we can take y arbitrarily large and make the dual problem strictly feasible. More specifically, taking z = e and y = γ(0, e), where γ is greater that the largest eigenvalue of Q + A * I (e), gives a positive definite matrix in (5).

Our semidefinite bounds for k-cluster

We now summarize some useful results about the nonlinear SDP bounds we use to solve k-cluster to optimality. These bounds are based on the general nonlinear SDP bounds of [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds for 0-1 quadratic problems leading to quasi-newton methods[END_REF], but are specialized here for the k-cluster problem.

For any matrix A, we denote by A F the Frobenius norm of A, which is defined as A F := A, A . For a real number a, we denote its nonnegative part by a + = max{a, 0}. We extend this definition to vectors and matrices as follows: for x ∈ R n , we define (x + ) i := (x i ) + , for i = 1, . . . , n, and for a symmetric matrix A, we define A + := U Diag(λ + )U T , where A has eigendecomposition A = U Diag(λ)U T , with eigenvalues λ ∈ R n and orthogonal matrix U ∈ R n×n . We denote similarly a -, x -, and A -.

Our first proposition is based on [MR13, Theorem 2]; we give the main ideas of the proof here to provide a convenient summary (using the notation of this paper) of the derivation of the semidefinite bounds used in this paper for the k-cluster problem.

Proposition 1. Let I be a set of triangle inequalities. For y ∈ R 2n+2 and z ∈ R |I| , let the matrix X I (y, z) ∈ S n+1 be defined by

X I (y, z) := [Q -B * (y) + A * I (z)] + . Let α > 0 and let F α I : R 2n+2 × R |I| → R be the function defined by F α I (y, z) := 1 2α X I (y, z) 2 F + b T y + e T z + α 2 (n + 1) 2 . ( 6
)
Then F α I (y, z) is a valid upper bound for problem (1); that is,

(KC) ≤ F α I (y, z), for all y ∈ R 2n+2 , z ∈ R |I| + . (7)
Proof 2. It is clear (see, e.g., [KMR12, Lemma 1]) that if X ∈ S n+1 satisfies diag(X) = e and X 0, then X F ≤ n + 1. Therefore, since α > 0, the following problem

(SDP α I ) maximize Q, X + α 2 (n + 1) 2 -X 2 F subject to B(X) = b, A I (X) ≥ -e, X 0, satisfies (SDP I ) ≤ (SDP α I ). ( 8
)
The Lagrangian of the (SDP α I ) problem is given by

L(X; y, z) := Q, X + α 2 (n + 1) 2 -X 2 F + y, b -B(X) + z, e + A I (X) = Q -B * (y) + A * I (z), X - α 2 X 2 F + b T y + e T z + α 2 (n + 1) 2 , where y ∈ R 2n+2 and z ∈ R |I| + . By [MR13,
Theorem 2] we have that F α I is the corresponding dual function; that is,

F α I (y, z) = max X 0 L(X; y, z).
Then, by weak duality we have

(SDP α I ) ≤ F α I (y, z), for all y ∈ R 2n+2 , z ∈ R |I| + . (9)
Combining inequalities (4), (8), and (9), we obtain inequality (7).

The bounds F α I (y, z) coincide, up to change of sign and notation, with the bounds Θ(λ, µ, α) of [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds for 0-1 quadratic problems leading to quasi-newton methods[END_REF]. For fixed α > 0 and set of triangle inequalities I, the best possible bound F α I (y, z) can be found by solving the problem minimize F α I (y, z) subject to y free, z ≥ 0.

The following proposition, based on [MR13, Theorem 2], gives us an important fact that supports the practical use of these bounds: the above problem is a convex and differentiable problem.

Proposition 2. The function F α I is convex and differentiable; its gradients are given by

∇ y F α I (y, z) = b - 1 α B (X I (y, z)) , and ∇ z F α I (y, z) = e + 1 α A I (X I (y, z)) . (10)
Proof 3. From the proof of Proposition 1, the function F α I can be interpreted as a dual function, which immediately implies its convexity. The differentiability of

F α I follows from [MR13, Theorem 2].
This smoothness allows us to use a quasi-Newton method that can handle bound constraints, such as L-BFGS-B [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF], to efficiently minimize F α I (y, z) over y and z ≥ 0. For the management of the different parameters, we closely follow the bounding procedure of [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] for Max-Cut: we minimize F α I (y, z) with increasing accuracy (the stopping tolerance ε of the quasi-Newton algorithm is driven to 0) while gradually adding inequalities and reducing α. The sketch of our bounding procedure is given in Algorithm 1; for more details, we refer to [KMR12, Algorithm 1].

Algorithm 1 Sketch of bounding procedure for k-cluster Input: Scalars α 1 > 0, ε 1 > 0, and initial set of triangle inequalities

I 1 = ∅ for k = 1, 2, . . . do Use L-BFGS-B to compute (y k , z k ) such that max b -B (X k ) ∞ , e + A I (X k ) -∞ < ε k , where X k ← 1 α k X I k (y k , z k ).
Remove inequalities I - k that are not active, and add some inequalities I + k that are violated by X k :

I k+1 ← (I k \ I - k ) ∪ I + k .
if the number of inequalities added

|I + k | is small then Decrease α k and ε k end if end for
We can extend the theoretical convergence analysis of [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] to obtain the following proposition, which corresponds, for the bounding procedure of k-cluster, to the Theorem 1 of [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] regarding the bounding procedure of Max-Cut.

Proposition 3. Let the sequence (α k , ε k , X k , y k , z k , I k ) k be generated by Algorithm 1. If (α k ) k and (ε k ) k both converge to zero, and ( X, ȳ, z, Ī) is an accumulation point of the sequence (X k , y k , z k , I k ) k , then the sequence of bounds converges to the classic semidefinite bound:

lim k→+∞ F α k I k (y k , z k ) = (SDPĪ). ( 11
)
The proof of Theorem 1 of [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] (and the lemmas and results used in the proof) extends easily to the k-cluster problem -except for Lemma 2 of [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF]. This lemma is in fact the obstacle to a straightforward generalization, since the proof given in [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] cannot be directly generalized here: it requires that the semidefinite relaxation is strictly feasible, which is not the case for the k-cluster problem, as discussed in the previous section. Let us derive here the convergence result without this argument.

Proof 4. Consider the subsequence that converges to ( X, ȳ, z, Ī). Since the I k are finite sets, the convergence of a subsequence to Ī means that there are infinitely many indexes k i of this subsequence such that I k i = Ī. By (6), we have F

α k i Ī (y k i , z k i ) = e T y k i + e T z k i + α k i ( X k i 2 + n 2 )/2.
Since α k → 0, the straightforward generalization for k-cluster of Lemma 3 of [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] implies that the sequence (F α k I k (y k , z k )) k converges. Passing to the limit in the bounds, we get

lim k→+∞ F α k I k (y k , z k ) = lim i→+∞ F α k i Ī (y k i , z k i ) = b T ȳ + e T z. ( 12
)
Let us now pass to the limit for the variables. We can repeat the arguments when passing to the limit in the proof of Lemma 2 of [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] (replacing diag and Diag by B and B * ): we get that X is feasible for (SDPĪ), that (ȳ, z) is feasible for (DSDPĪ), and that we have complementarity

Q -B * (ȳ) + A * Ī (z), X = 0 and zT (AĪ( X) + e) = 0. ( 13
)
Observe now that for all X feasible in (SDPĪ), we have

b T ȳ + e T z ≥ b T ȳ + e T z + Q -B * (ȳ) + A * Ī (z), X = Q, X + ȳT (b -B(X)) + zT (e + AĪ(X)) ≥ Q, X .
Moreover (13) gives that the above inequalities are equalities for X. This shows that the limit b T ȳ + e T z coincides with the optimal value of (SDPĪ). This proposition establishes the convergence and thus the consistency of the bounding procedure. More importantly, in practice the bounding procedure rapidly attains high quality bounds, which is the key property to be included in a branch-and-bound method. For illustration, we plot the convergence for our bounding procedure in Figure 1 on a k-cluster problem (n = 160, edge density 25%, k = 80). Such convergence curve is typical of what we have observed in our numerical experiments.

Branch-and-bound method for k-cluster

In this section, we describe our method for solving k-cluster to optimality. We list the main ingredients of our branch-and-bound implementation, and briefly compare them with the ones of [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF]. a) Bounds. We have described our bounding procedure in Algorithm 1 in the previous section. We now discuss the differences between the bounding procedure of Algorithm 1 and the one used in [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF]. With our notation, the bounds used in the branch-and-bound algorithm of [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF] are given by (see [MR12, Lemma 2]):

Θ(y, α) := α 2 [Q/α + B * (y)] + 2 F -αb T y+ α 2 (n+1) 2 , for α > 0. ( 14
)
Observe therefore that we have Θ(y, α) = F α ∅ (-αy) for all y and α > 0. In other words, the bounds used in [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF] essentially correspond (up to a scaling of the vector y) to the bounds that we use here, but with no triangle inequalities (I = ∅). The presence (and the efficient management) of the inequalities in Algorithm 1 is the key difference between our bounding procedure and the one of [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF]. Though it is well-known that triangular inequalities greatly improve the SDP bounds for k-cluster (see, e.g., [START_REF] Roupin | From linear to semidefinite programming: an algorithm to obtain semidefinite relaxations for bivalent quadratic problems[END_REF]), it is a computational challenge to incorporate them in the bounding in a tractable way; this has never be done for this problem in a context of exact resolution (neither by [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF] nor [START_REF] Pisinger | Upper bounds and exact algorithms for p-dispersion problems[END_REF][START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF]). In practice, Algorithm 1 gradually reduces the two tightness parameters α and ε to zero while incorparating blocks of promising inequalities. This is another difference with [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF] for which α and ε are kept fixed. Moreover, although both bounding procedures (ours and the one of [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF]) use quasi-Newton algorithms, here we have to use a quasi-Newton solver that can handle bound constraints, and we also have to run this solver several times during the computation of the bounds. b) Heuristics. We use the same two heuristic methods in the current branch-and-bound implementation as in the one of [START_REF] Pisinger | Upper bounds and exact algorithms for p-dispersion problems[END_REF] and [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF]: (a) For the initial feasible point, we use the classical two-step greedy heuristic for k-cluster, since it gives very good feasible solutions. (b) After running the bounding procedure on a subproblem having k ′ nodes added to the cluster, we add the remaining k -k ′ nodes having the largest fractional values.

c) Branching rules. In the current branch-and-bound implementation, we use the "difficult first" branching rule to decide which variables to branch on next: in the {0, 1} problem formulation, we select the variable having fractional value closest to 1 2 . This is the same branching rule that was used in [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF]. In addition, given a list of subproblems in the branchand-bound method, we must decide which subproblem to branch on next. We use the BOB Branch & Bound platform [START_REF] Cun | BOB: A unified platform for implementing branch-and-bound like algorithms[END_REF] that provides an easy and flexible way to implement a branch-and-bound algorithm. The BOB platform automatically handles the management of subproblems, and we select the subproblem with the weakest bound with the hope of making the most progress. In contrast, a basic depth-first traversal of the branchand-bound tree was used for the selection of subproblems in [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF].

Numerical results

As far as we are aware, the most challenging k-cluster test-problems publicly available are the ones used by [Bil05, BEP09, MR12], created by the graph generator rudy. The parameters of the instances of k-cluster are the size of the graph n, the value of k, and the graph density d. We take:

k = 1 4 n, 1 2 n, 3 4 n, d = 25%, 50%, 75%.
The instances are randomly generated as follows: given a density d ∈ [0, 1], a uniformly distributed random number ρ ∈ [0, 1] is generated for any pair of indexes i < j; if ρ > d, then w ij is set to 0, otherwise it is set to

• either to 1, generating an instance of pure k-densest subgraph problem,

• or to a integer in {0, . . . , 100} (or {-100, . . . , 100}), generating an instance of (generalized) heaviest k-subgraph problem.

We implemented our algorithm in C / FORTRAN and have used the Intel Math Kernel Library (MKL) for the eigenvalue computations. We report the computational results of our algorithm on the two above sub-families of k-cluster problems. In our tests we used a Dell Intel Xeon CPU E31270 3.40 GHz (using a single core) with 8Go of memory and running the Linux operating system.

Numerical results for k-densest subgraph

For the k-densest subgraph instances, we have n = 80, 100, 120, 140, 160, and five instances generated by the graph generator rudy for each set of parameters (n, d and k). We also have a set of random instances for n = 120, 140, 160 generated by Amélie Lambert and available online. Thus we have a total of 360 instances of the k-densest subgraph problem.

The instances with n = 80, 100 have already been used in previous articles about k-cluster, such as [Bil05, BEP09, MR12]; we call them the standard instances. We call the instances with n = 120, 140, 160 the larger instances. As far as we are aware, this paper is the first one where numerical results are reported for instances with n = 140 and n = 160. The instances with n = 120 were used by [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF], and their solver failed on some on them.

Our first experiment consists in comparing our method with the two best existing methods to solve the k-cluster problem to optimality:

1. the semidefinite method of [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF], which is a precursor of the method presented here using the same branch-and-bound structure but basic semidefinite bounds (without dynamic control of α, ε and without triangle inequalities), as described in the previous section; 2. the quadratic convex reformulation (QCR) of [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] which first convexifies the objective function of the problem (1) and then uses the state-of-the-art IBM/CPLEX mixed-integer convex quadratic solver (we use CPLEX 12.6 for the experiments). Since the best convexifying parameters of the initial problem are obtained by solving the semidefinite relaxation (SDP I ) (we use CSDP [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF] for the experiments), the QCR and our method use the same (tight) bound at the root node of the branch-and-bound tree. However, the tightness of the QCR bounds deteriorates when going down this the search tree.

We have implemented the three methods, tuned their parameters to reach the best performance, and run them on the same machine. Table 1 reports their comparison in terms of CPU time and number of nodes of the search tree to reach optimality for all the problems. Note that the reported times for QCR are the ones of CPLEX 12.6 (single thread). The dual variables needed to convexify the problem are obtained by CSDP [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF]; but the computing times of CSDP are not consider in the reported CPU time as they are usually unsignificant (about 0.4 seconds for n = 80 and 0.9 or n = 100).

The figures of Table 1 show that our algorithm clearly outperforms the two other methods. Regarding memory issues, the two first methods use a few amount of memory (less than 4 Mo) whereas the third featuring CPLEX uses up to 6 Go. In view of this first experiment on the standard instance, we know focus on our method and try to reach larger size.

For the second experiment, we consider the larger instances and report in Table 2 the average number of nodes and time required for our solver to solve each set of five problems. We emphasize that our solver does not need to visit a lot of nodes in the branch-and-bound search tree; for example, we have found that 55% of the problems with size n ≤ 120 are solved at the root of the branch-and-bound tree. Our algorithm is also able to solve unstructured k-cluster problems of sizes n = 140 and n = 160, for which, as far as we are aware, no numerical results have been reported in the literature. Finally we note that 96% of the 360 test problems in are solved within three hours.

There are significant differences on the performance of the algorithm for problems of the same size. In particular, we notice that the number of nodes and total time tends to increase when d increases from 25%-75%. This is due to the presence of many near-optimal k-clusters (and more generally many kclusters can have very similar values) for larger density graphs. Consequently, even if the bound is tight, it is hard to prune nodes in the search tree since the evaluations of many nodes are almost the same. The bottomline is that n = 160 is the largest size of unstructured k-cluster problems to be solved to optimality within a reasonable amount of time on a single-threaded machine.

Numerical results for heaviest k-subgraph

For the (generalized) heaviest k-subgraph instances instances, we have n = 120, 140, 160 and five instances generated by the graph generator rudy for each set of parameters (n, d, and k) and for positive weights in 0, . . . , 100 and for integer weights in -100, . . . , 100. This makes a total of 270 instances. As above, we report the average number of nodes and time required to solve each set of five problems in Table 3.

The computing time and number of nodes reported in the two Tables 2 and 3 are comparable. The main difference between the two tables is for dense graphs with large k (look for example at the last lines of the two tables). For dense graphs with large k, unweighted instances have several optimal solutions (thus, it is easier to find an optimal solution). When k is smaller, the instances in the weighted and unweighted cases are more difficult and comparable with each other (look, for example, at the first line with n = 160, in the two tables).

The bottonline is that the nature of the weights (positive, negative, or 0-1) do not change the order of magnitude of computing times and the number of nodes. This confirms the folklore knowledge that the instances with {0, 1} weights (k-densest subgraph instances) are as hard to solve as the weighted ones (heaviest k-subgraph instances). This is probably the reason why [Bil05, BEP09, MR12] only considered these instances. Such numerical observations can be related to theoretical results in complexity theory for similar optimization problems in graphs (see, e.g., [START_REF] Crescenzi | On weighted vs unweighted versions of combinatorial optimization problems[END_REF]).

Online complementary material and on-going developments

We have made the following complementary material available online:

• the entire dataset of problems,

• the full numerical results of our tests,

• a web interface for the solver.

To access this material, we invite the interested reader to visit the BiqCrunch website: http://lipn.univ-paris13.fr/BiqCrunch. and linear constraints. Since the theory has already been presented for the general case [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds for 0-1 quadratic problems leading to quasi-newton methods[END_REF], our work essentially consists of modeling and computer implementations. For example, during the review process of this paper, we have done numerical experiments with the max-clique problem (also called "max independent set problem"), for which there exist many publicly available instances (DIMACS instances). These complementary results are and will be be posted on the web site. We want to point out that it is for the k-cluster problem that we have got the most spectacular computational results, as reported in this paper. This comes probably from the fact that k-cluster is harder than other problems -in particular, harder than max-clique, both in the sense of computational resolution and in the sense of complexity theory (max-clique reduces to kcluster in bipartite graphs [START_REF] Corneil | Clustering and domination in perfect graphs[END_REF]). Thus we do not present other results in this paper and instead we refer the interested reader to the website for further research and developments.

Conclusions

In this paper, we have presented an improved branch-and-bound algorithm to solve k-cluster problems to optimality. Our method is based on the previous work of [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF] and [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] and here we have summarized the main ideas from these papers and have highlighted the differences. The two main reasons why the method presented here is much better than the results reported in [START_REF] Malick | Solving k-cluster problems to optimality with semidefinite programming[END_REF] are the introduction of triangle inequalities and the strategy of reducing α and ε in the semidefinite bounding procedure presented in Algorithm 1.

The main contribution of this paper is the extensive numerical experiments on benchmark k-cluster problems, including large-scale problems that are beyond the ability of previous exact methods to solve within a reasonable amount of time. We were able to solve for the first time a set of hard k-cluster instances of size n = 160, advancing our ability to solve instances from the previous limit of n = 120.

Our other contributions are: adapting the branch-and-bound method and bounding procedure of [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] to solve k-cluster problems; extending the convergence results of [START_REF] Krislock | Improved semidefinite bounding procedure for solving max-cut problems to optimality[END_REF] in spite of the lack of strict-feasibility of the SDP relaxation of Lemma 1; comparing practical complexity for different subfamilies the of k-cluster problem; making the benchmark dataset of kcluster problems and entire numerical results available online together with a web interface for our solver.

Since the semidefinite relaxation (2) of the k-cluster problem is not strictly feasible, we will also consider in our future work the use of semidefinite facial reduction (see, e.g., [START_REF] Ramana | Strong duality for semidefinite programming[END_REF]) with the hope of further improving the time to solve k-cluster problems to optimality.

Figure 1 :

 1 Figure 1: Time/bound plot of the Algorithm 1 bounding procedure on problem PB3 n160 d025 k80 which has optimal value 1026.

Table 1 :

 1 For the standard k-cluster problems: for the three best existing methods, we compare the number of nodes and CPU time (s), averaged over five instances for each triple (n, k, d).

		(b) our method	(c) [MR12]	(d) [BEP09]
	n k d(%)	nodes time	nodes time	nodes time
	80 20 25	3.4 3.2	11650	94.5	170658	39.9
	50	7.4 6.1	41857 323.1	536648 125.0
	75	13.8 9.5	102948 1002.4	1827452 407.4
	40 25	1.4 1.1	1544	13.1	26597	8.6
	50	1.0 0.8	2806	24.6	34148	9.6
	75	6.6 8.0	19789 195.7	231620	55.2
	60 25	1.0 1.1	148	1.3	946	0.6
	50	1.0 0.8	302	2.7	5128	3.1
	75	1.0 0.7	1123	11.2	5754	3.4
	mean	4.1 3.5	20241 185.4	315439	72.5
	100 25 25	20.6 31.3	127901 2207.1	5680415 1882.8
	50	35.0 41.6	303648 5543.0	19164583 6684.9
	75	30.6 32.3	1180710 19661.2	44336562 14275.8
	50 25	3.4 5.3	9328 164.9	415340 153.6
	50	25.4 46.1	211308 3923.5	5156390 2182.2
	75	1.4 2.1	27099 514.2	514822 203.1
	75 25	1.0 1.9	455	8.0	10261	5.1
	50	1.8 3.8	2018	39.7	108962	37.3
	75	1.0 1.6	1958	38.7	14956	6.5
	mean	13.4 18.5	207158 3566.7	8378032 2825.7

Table 2 :

 2 For the larger problems: the number of nodes and CPU time, averaged over five instances for each triple (n, k, d). Entire set of results available online at http://lipn.univ-paris13.fr/BiqCrunch/results.

		(b) rudy instances	(c) other instances
	n k d(%)	nodes time (s)	nodes time (s)
	120 30 25	119.4	315.8	16.6	36.0
	50	194.2	425.1	39.4	83.1
	75	422.2	889.8	62.2	119.5
	60 25	59.8	198.5	12.2	28.3
	50	85.8	263.2	27.4	69.4
	75	43.0	143.0	41.8	93.0
	90 25	1.8	6.8	1.0	2.9
	50	22.2	96.9	1.0	2.7
	75	1.0	3.0	1.0	2.6
	140 35 25	366.2	1165.8	131.0	445.1
	50	1063.4	2888.6	383.0	964.7
	75	1558.6	4079.5	485.0	1279.8
	70 25	134.2	543.0	54.2	216.5
	50	780.6	3035.3	298.6	1133.0
	75	52.2	202.9	155.8	571.7
	105 25	2.6	14.1	1.4	7.5
	50	11.0	61.5	7.4	39.8
	75	6.6	34.8	3.0	19.1
	160 40 25	744.6	2856.4	235.4	1023.6
	50	11325.4 37565.2	858.6	3280.1
	75	8050.6 26302.6	1132.2	3997.8
	80 25	395.4	1835.2	73.4	401.6
	50	993.4	4654.5	479.8	1908.6
	75	3829.0 18653.9	1425.0	6288.9
	120 25	31.4	219.8	1.4	9.5
	50	17.4	143.4	2.2	16.7
	75	9.8	82.2	4.2	31.7

Table 3 :

 3 For graphs with weights: the number of nodes and CPU time, averaged over five instances for each triple (n, k, d). Entire set of results available online at http://lipn.univ-paris13.fr/BiqCrunch/results.

Acknowledgment. We thank the referee for their pertinent comments and suggestions which helped us improve our numerical comparisons, as well as the convergence analysis. We acknowledge the support of the grant "ANR GeoLMI" and of the CNRS "Mastodons" grant "Gargantua".


Laurent Boudin

Francesco Salvarani

Opinion dynamics: kinetic modelling with mass media, application to the Scottish independence referendum

Introduction

Forecasting the opinion evolution with respect to a binary question is crucial in many situations. A typical example consists in the anticipation of a referendum results of competition by using poll data from surveys held a few weeks before the vote. Our work lies in the sociophysics domain, which was born in the 80s. The reader can check [START_REF] Galam | Sociophysics: a physicist's modeling of psycho-political phenomena (understanding complex systems[END_REF] and the references therein to know more about the topic. In particular, one can find there some discussions about the French referendum on the European constitution in 2005.

In this article, we give a contribution to this problem by studying a mathematical model based on a kinetic approach. We perform the following hypotheses. Assumption 1. The number of individuals in the population is constant. Of course, this is only relevant for short-term forecasts, as during a referendum campaign. Assumption 2. The community is fully interconnected and all the members can somehow discuss with each other. By tuning the cross section of the model, however, it is possible to generalize our model to situations where the interaction probability between individuals depends on their own convictions.

Assumption 3. There exists a number of mass media which can have an effect on the individuals. This assumption is natural in most countries, since media are both observers and opinion carriers.

Let us more focus on the mass media, whose involvement is the most significant contribution of this work. Although the press freedom appears as a safeguard in a democracy [START_REF] Lichtenberg | Democracy and the mass media: a collection of essays[END_REF], the media can also be used to influence the public opinion, by taking advantage of their possibility of easily reaching a wide audience and favour some interests, ideas or arguments inside a population [START_REF] Coxall | Human manipulation, a handbook[END_REF]. In order to influence the population, many different tactics have been developed. Among them, we quote the distraction, in different forms (which is based on the assumption that the public has a limited attention span), the appeal to consensus and the fear mongering [START_REF] Ellul | Propaganda: the formation of men's attitudes[END_REF]. These techniques are used either separately or combined between them, and could be applied to many different contexts. However, they need some common factors in order to be successful.

Mass media use diversified media technologies in order to reach a large audience and transfer concepts, ideas, etc. Among them, we can cite broadcast media (radio, television), print media (newspapers, books), outdoor media (posters), digital media (Internet, mobile mass communication).

However, mass media often simply spread out the opinions which they support into the population, without any other interest than the maximisation of their audience.

These two different strategies cannot be modelled in the same way. Manipulating media can be described as entities with a given external opinion, which aim to spread their own opinion inside the population. Wide audience media should be described with a more intricate approach: the opinions they carry can depend on the opinion distribution inside the population. Hence, in this latter case, the influence of the population on the media opinion has to be part of the model itself.

In this article, we only consider media with an a priori given opinion. Note that this case can be seen as a first approximation of the wide audience media, where the media opinion may not remain constant, but does not depend on the opinion distribution.

We here use a kinetic approach to describe and forecast the evolution of a system under the aforementioned effects. In this framework, the information on the system is contained in a distribution function whose time evolution is governed by a partial differential equation with integral operators. This strategy is based on sophisticated mathematical tools and its interest is apparent when the number of individuals is very large, since it allows to handle a collective behaviour.

Kinetic equations have been used to model social phenomena since the early 90s, when Helbing studied behavioural changes by using Boltzmannlike equations [START_REF] Helbing | Stochastic and Boltzmann-like models for behavioral changes, and their relation to game theory[END_REF][START_REF] Helbing | Boltzmann-like and Boltzmann-Fokker-Planck equations as a foundation of behavioral models[END_REF][START_REF] Helbing | A mathematical model for the behavior of individuals in a social field[END_REF]. Subsequently, this approach has been the basis of several works, see the review article [START_REF] Boudin | Modelling opinion formation by means of kinetic equations[END_REF] and the references therein. Note that, in most part of the literature based on the kinetic approach, the key phenomenon is compromise (see, for example, [START_REF] Ben-Naim | Bifurcation and patterns in compromise processes[END_REF][START_REF] Ben-Naim | Unity and discord in opinion dynamics[END_REF][START_REF] Toscani | Kinetic models of opinion formation[END_REF][START_REF] Düring | Boltzmann and Fokker-Planck equations modelling opinion formation in the presence of strong leaders[END_REF]).

The model we here investigate also owns this binary interaction feature, by supposing that individuals follow the rule proposed by the authors in [START_REF] Boudin | A kinetic approach to the study of opinion formation[END_REF]. By taking account the influence of mass media on the population, we are able to improve our previous results in [START_REF] Boudin | A kinetic approach to the study of opinion formation[END_REF] and recover some realistic behaviours.

The opinion variable of our model runs in a continuous way between two extreme values, from -1 to +1. The model describes the time evolution of a distribution function f (depending on both opinion and time), which represents the density of individuals having an opinion about the binary question. It reflects the opinion formation process, which is not necessarily of binary type, since the doubt and the partial agreement are often predominant feelings.

The media action on the system is modelled by a kinetic operator which has a structure similar to the media operator introduced in [START_REF] Boudin | Kinetic model for multidimensional opinion formation[END_REF], but with a different nonlinear effect in the post-interaction opinion. The structure of this linear term (where the linearity is intended only with respect to f ) is well adapted to the coupling with the self-thinking.

At the end of the opinion formation process, the opinions must be translated into a decision. This second step is not tackled here, since it is often the result (at least in the case of reasonable and rational individuals) of a game-theoretical approach, especially when those individuals have intermediate opinions.

Since we are working with a continuous opinion variable, it is difficult, in general, to get comparisons with real data. Indeed, polls usually have a binary or ternary structure, since they are interested in forecasting the result of the final choice with respect to the binary question of the referendum. However, some more structured surveys have been built on a more complete scale (usually from 1 to 10). These polls can be a good tool for comparisons, which are essential to somehow validate our model. The last part of our study is precisely devoted to the qualitative comparison between the results of our model and three surveys performed by the polling corporation ICM Unlimited1 about the Scottish independence referendum, that took place in Scotland on September 18, 2014.

The article is organized as follows: in the next section, we describe the model and its mathematical structure. The existence of solutions is investigated in Section 3, whereas the numerical simulations are performed and commented in Section 4. Eventually, Section 5 focuses on the Scottish independence referendum: we relate our numerical results to some existing polls obtained before the referendum itself. It provides some kind of validation to our model, and already suggests possible improvements.

Kinetic model

In this section, we briefly describe the model we investigate here. As already stated, it appears as an extension to the one previously introduced in [START_REF] Boudin | A kinetic approach to the study of opinion formation[END_REF]. We refer to that article for a full description of the bilinear integral operator defining the binary interaction between the agents.

In what remains, Ω denotes the open interval (-1, 1). The variables of the model are opinion x, a continuous variable belonging to Ω, and time t ∈ R + . The opinion variable x describes the degree of agreement with respect to a binary question (for example, a referendum). In particular, x = ±1 are the two extreme answers to the question, i.e. "yes" or "no" without reserve, whereas any intermediate value between -1 and +1, 0 excluded, means that the corresponding agent partially agrees with the opinion labelled with the same sign, with a degree of conviction which is proportional to |x|. The value x = 0 means that the agent is undecided.

The community is described by means of the distribution function f := f (t, x), defined on R + × Ω, whose time evolution is governed by a kinetic equation. The main ingredients of the kinetic model act on two different levels. The first one is the description of the microscopic active phenomena, in this case the self-thinking, the binary opinion exchange between the individuals of the population and the effect of mass media on the single individual. The second level governs the time evolution of the distribution function, which is induced by the operators which take the microscopic phenomena into account. Since the model is of kinetic type, we borrow the language of kinetic theory. For instance, the term collision means an interaction with exchange of opinions, that gives, as a result, a modification of the agents opinions.

2.1. Self-thinking. The self-thinking phenomenon is described by a diffusion operator obeying to a non-uniform Fourier law, with Fourier term α = α(x). This term quantifies the possibility that people may change their opinion through personal reasoning. In particular, we assume that the Fourier term satisfies the following properties.

Definition 2.1. Let α : Ω → R be a nonnegative function of class C 1 ( Ω). We say that α is admissible if • α(x) = α(-x) for any x ∈ Ω, • α(-1) = α(+1) = 0.
2.2. Binary opinion exchanges. Let x, x * ∈ Ω the opinions of two individuals of the population before interacting, and x ′ , x ′ * ∈ Ω the corresponding opinion after the binary exchange.

We suppose that the individuals of the population are of conciliatory type [START_REF] Boudin | A kinetic approach to the study of opinion formation[END_REF][START_REF] Boudin | Conciliatory and contradictory dynamics in opinion formation[END_REF]: it means that they have a natural trend to reach a compromise with each other. From a quantitative point of view, the opinions after interaction are modified in order to get closer to the average opinion before the interaction and, at the same time, stronger opinions are less attracted towards the average than weaker ones. More precisely, the collision rule adopted in the article is the following:

(1)

x ′ = x + x * 2 + η(x) x -x * 2 , x ′ * = x + x * 2 + η(x * ) x * -x 2 .
The function η : Ω → R is called the attraction function, and satisfies the following properties.

Definition 2.2. A function η ∈ C 1 ( Ω) is admissible if 0 ≤ η < 1, η ′ (x) > 0
for all x > 0, and the Jacobian J(x, x * ) of the collision mechanism (1) is lower bounded by a positive constant, i.e. there exists J 0 > 0 such that J(x, x * ) ≥ J 0 for any x, x * ∈ Ω.

The collision rule (1) governs the interaction between individuals and the corresponding exchange of opinions, through a collisional integral of Boltzmann type.

2.3.

Interaction with mass media. In the model, we take into account the existence of m mass media, m ≥ 1. They can influence the population by sharing their opinion about the referendum. The effects of the media are modelled thanks to an interaction with a given background. To define the characteristics of each media M j , 1 ≤ j ≤ m, we need three quantities: its time-depending strength θ j : R + → R + , its time-depending opinion X j : R + → Ω, and its attractiveness q j : Ω2 → R.

The individuals of the population are still considered as conciliatory regarding the media. The evolution of an individual's opinion follows a rule, similar to the binary collision mechanism defined in the previous subsection, but with the difference that the opinion of the media is not influenced by the opinion of the population. The main idea of the collision mechanism consists in supposing that the opinions after the interaction will be modified in order to approach the average opinion before the interaction and, again, that stronger opinions are less attracted towards the media opinion than weaker ones. Hence, if we denote ny xj the post-collisional opinion of an individual of pre-collisional opinion x, for each mass media M j , j ∈ [1, . . . , m], we write

(2) xj = x + X j (t) 2 + q j (x, X j (t)) x -X j (t) 2 ,
where q j satisfies the following properties.

Definition 2.3. For any j, the media attractiveness

q j ∈ C 0 ( Ω × Ω) is admissible if 0 ≤ q j ≤ 1 for any x ∈ Ω.
Moreover, it seems reasonable, from the modelling point of view, to impose that the media attractiveness equals 1 when |x| ≃ 1: it means that individuals with a strong opinion are less influenced by the media effects than agents with a weaker opinion.

2.4. Time evolution of f . The model governing the evolution of the system is a partial differential integral equation, which can be more naturally written in a weak form. Indeed, the existence of a pre-collisional opinion pair generating a post-collisional opinion pair (x, x * ) through the collision rule (1) is not guaranteed, in general. The model hence reads as follows. Let T > 0 and f in ∈ L 1 (Ω) a nonnegative function. The unknown f satisfies, for almost every t ∈ [0, T ] and in a distributional sense in t,

(3) d dt Ω f ϕ dx = Ω ∂ x (α∂ x ϕ) f dx + Q(f, f ), ϕ + m j=1 L j (f ), ϕ for any test-function ϕ ∈ C 2 ( Ω), with initial condition (4) f (0, x) = f in (x), x ∈ Ω.
Operator Q translates the effects of the binary interactions between the agents, and can be written in different forms. A key ingredient of this bilinear (with respect to f ) operator is the cross-section. This nonnegative quantity, denoted by β, governs the probability that an exchange of opinions can occur. We here assume that β is a positive constant. This is the simplest possible assumption. It implies that the probability of interaction between two individuals depends neither on their respective opinions, nor on time, and that there are no forbidden exchanges.

The weak form used in this work is

(5) Q(f, f ), ϕ = β Ω 2 f (t, x) f (t, x * ) ϕ(x ′ ) -ϕ(x) dx * dx.
Note that the collision rule (1) then only appears in the arguments of the test-functions. Operators L j , 1 ≤ j ≤ m, describe the media effect on the population, and they have the following form:

(6) L j f, ϕ = θ j (t) Ω f (x) (ϕ(x j ) -ϕ(x)) dx.
Note that those operators linearly depend on f .

Main mathematical properties

It is easy to deduce some basic mathematical properties, such as the conservation of the number of individuals in the population. Proposition 3.1. Let f = f (t, x) be a nonnegative weak solution of (3)-( 4), with a nonnegative initial datum f in ∈ L 1 (Ω). Then we have

f (t, •) L 1 (Ω) = f in L 1 (Ω)
for a.e. t ≥ 0.

Proof. This is a straightforward consequence of (3)-( 6), with test-function ϕ ≡ 1. Some mathematical properties of the equation are a consequence of the structure of the bilinear term. The following result holds, the proof of which is given in [START_REF] Boudin | A kinetic approach to the study of opinion formation[END_REF].

Lemma 3.2. Let f (t, •) ∈ L 1 (Ω). Then both Q + (f, f )(t, •) and Q(f, f )(t, •)
belong to L 1 (Ω), and we have, for a.e. t > 0,

Q + (f, f )(t, •) L 1 (Ω) ≤ 2β 1 -max η f (t, •) 2 L 1 (Ω) , (7)
Q(f, f )(t, •) L 1 (Ω) ≤ 2 1 -max η + 1 β f (t, •) 2 L 1 (Ω) . (8)
Existence of a nonnegative solution to (3)-( 4) can be deduced by construction, starting from the following result, again in [START_REF] Boudin | A kinetic approach to the study of opinion formation[END_REF].

Proposition 3.3. Consider the initial-boundary value problem for the unknown

v = v(t, x), x ∈ Ω and t ∈ [0, T ], (9) v t -[α(x)v x ] x + µv = g, µ ≥ 0, with initial condition (10) v(0, •) = v in
and boundary conditions

(11) lim x→±1 α(x)v x (t, x) = 0 a.e. t,
where v in ∈ L 1 (Ω), g ∈ C([0, T ]; L 1 (Ω)) are nonnegative functions. Then (9)-( 11) admits a unique solution v ∈ C 0 ([0, T ]; L 1 (Ω)), and v is nonnegative.

The procedure to build a nonnegative weak solution is based on the monotonicity properties of our problem. Theorem 3.4. Let f in a nonnegative function in L 1 (Ω). Then there exists a nonnegative weak solution f ∈ L ∞ (0, T ; L 1 (Ω)) to (3)-(4).

Proof. Set ρ = Ω f in (x * )dx * .
Consider the sequence (f n ) n∈N of functions, defined as solutions to ( 12)

Ω ∂ t f n+1 ϕ dx - Ω ∂ x (α∂ x ϕ) f n+1 dx +   βρ + j θ j (t)   Ω f n+1 ϕ dx = Q + (f n , f n ), ϕ + j θ j (t) Ω f n ϕ(x j ) dx,
with f 0 ≡ 0, for all ϕ ∈ C 2 ( Ω), satisfying the initial and boundary conditions

f n (0, •) = f in , n ≥ 1, lim x→±1 α(x)f n x (t, x) = 0 for a.e. t, n ≥ 1.
Thanks to Lemma 3.2, we can apply Proposition 3.3 and deduce, by induction, that f n exists, belongs to C 0 ([0, T ]; L 1 (Ω)) and is nonnegative. Then, choosing ϕ = 1 in ( 12), we obtain

d dt Ω f n+1 dx +   βρ + j θ j (t)   Ω f n+1 dx = β Ω f n dx 2 + j θ j (t) Ω f n dx.
Therefore, by finite induction, we immediately get

Ω f n dx ≤ ρ, n ≥ 1.
Moreover, by applying the same strategy as in [START_REF] Boudin | A kinetic approach to the study of opinion formation[END_REF], we can prove that (f n ) is a non-decreasing sequence.

By monotone convergence, we finally deduce the existence of f as the limit of (f n ) in L ∞ (0, T ; L 1 (Ω)). In order to check that f satisfies (3), we write

T 0 Ω f n+1 ϕ(x) ψ(t) dx dt - T 0 Ω α(x) ϕ ′ (x) ′ f n+1 ψ(t) dx dt +   βρ + j θ j (t)   T 0 Ω f n+1 ϕ(x) ψ(t) dx dt = β T 0 Ω 2 f n (t, x)f n (t, x * ) ϕ(x ′ ) ψ(t) dx dx * dt + j T 0 Ω θ j (t) f n ϕ(x j ) ψ(t) dx dt.
The only small difficulty lies in the fourth integral. It also converges since

Ω 2 |f n (x)f n (x * )-f (x)f (x * )| |ϕ(x ′ )|dxdx * ≤ 2ρ f n (t, •)-f (t, •) L 1 ϕ L ∞ .
That ends the proof of Theorem 3.4.

Numerical tests

The architecture of the numerical method used in the simulations treats separately the kinetic part and the diffusion term by means of a splitting technique. More details about the numerical scheme are given in [START_REF] Boudin | A kinetic approach to the study of opinion formation[END_REF], except for the media numerical treatment. The mass media are a set of invariant particles which carry the media opinion and interact with the distribution function f , chosen by a sample over the whole set of macro-particles, by obeying to the collision rules [START_REF] Ben-Naim | Unity and discord in opinion dynamics[END_REF].

In all the computations, the Fourier term is α(x) = 0.05(1 -x 2 ) 1/3 , which is in agreement with the hypotheses of Definition 2.1, the collision frequency β is set to 50, and the attraction function is given by η(x) = 0.25(1 + x 2 ), which satisfies the constraints of Definition 2.2.

The media attractivenesses q j are crucial. In our tests, we shall consider two different types of media influence: a global media influence function and a media influence which equals 1 outside an interval [-δ, δ], 0 < δ < 1. This last choice translates the idea that strong opinions are not affected by the influence of the media, whereas weaker opinions are sensitive to the effect of mass media.

In what follows, we investigate some relevant test cases, which give an overview of the quantitative features of the model. First of all, we study the influence on a population of a unique media with a fixed opinion, and two different initial conditions: a balanced initial datum and an unbalanced one. Subsequently, we discuss the interactions between two groups of balanced media, again with fixed opinions. This last result is compared in Section 5 to the opinion dynamics for the Scottish independence referendum.

4.1.

Effect of a unique media on the opinion formation. The first example studies the behaviour of the initial value problem (3)-( 4), with m = 1. In all tests of this subsection, we suppose that θ 1 is constant equal to 50.

The media attractiveness is also constant: q 1 = 0.5. In Figure 1, we can check that the asymptotic states of the system with media opinion X 1 (t) = -0.9, for both initial data f in 1 = 0.5 (plot on the left) and f in 2 = 1 [0,1] (plot on the right).

We note that the asymptotic states are the same, and of course do not depend on initial data. They are driven by the media opinion, even if the population has an initial opinion with opposite sign with respect to the media opinion. Moreover, both situations show that self-thinking prevents the system from a complete adherence to the media opinion.

In Figure 2, we can see the sensitivity of the model with respect to the media attractiveness q 1 . We here choose a regularized version of

q 1 = 1 - 1 2 1 [-0.5,0.5] .
Note that, since we are working on a discrete grid, we can numerically identify the discretization of q 1 and the discretization of its regularized version.

This media attractiveness means that only individuals with a weak opinion, in the range [-0.5, 0.5], are influenced by the media opinion.

In order to compare the results with the previous case, we choose the same set of initial datum and media opinion as in the second example, i.e. f in 2 = 1 [0,1] and X 1 (t) = -0.9. After a transient period of time, the media has, again, enough strength to drive the whole population towards opinions whose sign agrees with its opinion, even if 50% of the population cannot initially be influenced by the media opinion.

This behaviour is due to the binary interactions between the agents. Since the whole population is composed of conciliatory individuals, the effect of the binary interactions leads to compromise. Hence, the fraction of the population which keeps an opinion outside (0.5, 1] decreases towards 0. On the other hand, since the agents with opinions inside the interval [-1, -0.5) are not influenced by the media, the peak of the distribution is not centred at the media opinion any more, and is closer to 0.

4.2.

Interaction between two media with fixed opinions. We conclude by pointing out how a plurality of media is a factor that allows pluralism of opinions and the splitting of the population in two well-defined groups.

In order to recover this behaviour, we now choose m = 2, θ 1 = θ 2 = 50, and X 1 (t) = -X 2 (t) = 0.9.

In Figure 3, we plot the asymptotic state of the distribution function, with initial datum f in 1 = 0.5 and both media attractivenesses given by q j = 1 if |x -X j | > 1 and q j = 0.5 otherwise, for j = 1, 2. After a transient period, we note the formation of two peaks centred at each media opinion ±0.9. The fraction of agents with a weak opinion is very low. It is a straightforward consequence of the fact that all individuals of the population are under the media influence.

In Figure 4, we plot the asymptotic state of the distribution function, with the same initial datum as in Figure 3, namely f in 1 = 1/2. The media attractivenesses are given by q j = 1 if |x -X j | > 0.5 and q j = 0.5 otherwise, for j = 1, 2.

In this case, the individuals with opinion close to 0 are not influenced by any media: they can change their opinions only through self-thinking and binary interactions.

We note the formation of two strong sub-groups with asserted opinions, but also the emergence of a fraction of undecided persons (who may have a higher propensity to abstention). This behaviour is typical of some referendum campaigns, such as the case discussed in the next section.

Comparison with polls: the Scottish case

On September 18, 2014, Scotland was called to answer the independence referendum question: "Should Scotland be an independent country?". The turnout of 84.6% was the highest one recorded for a referendum in the United Kingdom since the introduction of universal suffrage. This datum is enough, by itself, to prove the importance that Scotland attributed to the The referendum was preceded by hundreds of surveys, whose aim was the prediction of the final result. Those surveys took into consideration many different aspects related to the independence referendum question. Among them, we have considered three surveys for the validation phase, which have the advantage to offer multiple answers and allow the comparison with our model, based on a continuum of opinions. The precise data of these polls are the following:

(1) July The question was identical in the three polls: "Can you please say where you are on this scale regarding Scotland becoming independent?". The answers were modulated on a scale, between 1 and 10, where

• 1 means "completely against an independent Scotland",

• 5 means "neither for nor against" (positioned between 5 and 6),

• 10 means "completely for an independent Scotland".

The results of the three surveys are all plotted in Figure 5. As the referendum day came, there was a clear bipolarization of the population and the emergence of two opposite fields with approximately the same magnitude. The number of undecided and agents with mild opinion decreased, and the media widely spread the arguments of both opposite parties.

The qualitative agreement between the surveys history and the results in Figure 4 is quite good: both peaks are centred at the extreme values (or very close to them). We recognize moreover the existence of some undecided people, which have been reported in the polls (there is a local maximum around the value 5).

The main differences between our simulations and the three surveys is given by the behaviour of moderately decided individuals (located in states 2 to 4 and 7 to 9). Indeed, the polls report that a small, but not negligible, fraction of polled people is located in these intermediary states, whereas, in our numerical simulation, the fraction of the population far from the peaks is very small. However, it is difficult to draw some conclusions from this consideration since the magnitude of the population in each of the intermediary states is often below 5%, and therefore it is quite close to the survey error margin.

Conclusion

We presented a kinetic model to describe an interconnected population which must choose an option in the framework of a binary question. The agents can modify their opinion by means of spontaneous self-thinking, discussions with each other, and influence of mass media. The model has been investigated both mathematically and numerically. The comparison between the model and the results of three surveys on the Scottish independence referendum of 2014 provides a good qualitative agreement.

Figure 1 .

 1 Figure 1. Effect of a unique media for two different initial data.

Figure 2 .

 2 Figure 2. Effect of a unique media with finite influence range.

Figure 3 .

 3 Figure 3. Pluralism of media and pluralism of opinions: global media influence

Figure

  Figure 4. Pluralism of media and pluralism of opinions: local media influence.

Figure 5 .

 5 Figure 5. Results of the three ICM surveys.

L

  .B.: Sorbonne Universités, UPMC Univ Paris 06, UMR 7598 LJLL, Paris, F-75005, France & CNRS, UMR 7598 LJLL, Paris, F-75005, France & Inria, Équipe-projet Reo, BP 105, F-78153 Le Chesnay Cedex, France E-mail address: laurent.boudin@upmc.fr F.S.: Dipartimento di Matematica, Università degli Studi di Pavia, Via Ferrata 1, I-27100 Pavia, Italy E-mail address: francesco.salvarani@unipv.it

  4. Pluralism of media and pluralism of opinions: local media influence. referendum question. The official results, certified by the Chief Counting Officer, from the Electoral Management Board for Scotland (EMB), are the following:

	Yes	1,617,989	44.7%
	No	2,001,926	55.3%
	Valid votes	3,619,915 99.91%
	Invalid or blank votes	3,429	0.09%
	Total votes	3,623,344 100.00%

  11, 2014: ICM survey for Scotland on Sunday with 1002 respondents, aged 16 and older; (2) August 13, 2014: ICM survey for Scotland on Sunday with 1005 respondents, aged 16 and older; (3) September 16, 2014: ICM survey for The Scotsman with 1175 respondents, aged 16 and older.

See http://www.icmunlimited.com/media-centre/polls/

This work was partially funded by the 2013-2014 French CNRS PEPS HuMaIn project Kindymo, headed by L. Boudin, the French ANR-13-BS01-0004 project, Kibord, headed by L. Desvillettes, and the Italian national institute of higher mathematics (INDAM), GNFM project Kinetic description of collective phenomena: theory and numerical simulations, headed by F. Salvarani.