Chapter 1

Interactive Scores: Past, Present and Future

Concurrent Constraints Time Model

Process calculi such as the Non-deterministic Timed Concurrent Constraint (ntcc) [NIE 02] calculus has been used in a variety of musical applications (see Chapter ??). An advantage of ntcc is that process synchronization is achieved by adding or deducing constraints from a constraint store, thus is declarative.

In this section we present a model of interactive scores based upon ntcc. The model is inspired on a previous model based upon ntcc that do not consider hierarchy nor the fact that execution can continue even if interactive events are not launched [ALL 06]. The model is also close is spirit to the petri nets model. Regrettably, in petri nets is difficult to model global constraints such as the maximum number of simultaneous temporal objects and temporal reductions during execution.

Another related model is Tempo, a formalism to define declaratively partial orders of musical and audio processes [RAM 06]. Unfortunately, Tempo does not allow us to express choice when multiple conditions hold, simultaneity nor to perform an action if a condition cannot be deduced.

Ntcc model

We use ntcc to the express operational semantics of interactive scores. In order to define operational semantics, we need to transform an interactive score into a graph Chapter written by Myriam DESAINTE-CATHERINE and Antoine ALLOMBERT and Mauricio TORO.

where the vertices are the start and end points of temporal objects and arcs represent the delays among them. We define the graph as g = (V, A, lV, lA) where V is the set of vertices, A the set of arcs, lV a function that assigns labels to vertices, lA a function that assigns labels to arcs. A vertex is labeled with the type of point (a start point, end point or interactive point) and the temporal object associated to it. An arc is labeled with its duration. The graph is reduced by removing zero-delay arcs and representing several points in the same vertex.

Absence of zero delays simplifies the definition of operational semantics because we do not have to synchronize two processes to occur at the same time. To remove a zero-delay arc between a vertex a and b, we delete a and we connect all its successors and predecessors to b. We also combine the label of a with the label of b, this means that a vertex may represent the start or end of several points.

Points

We have two type of points: interactive (iP oint) and static points (sP oint). Process Launch i updates the variable launched i with true and persistently assigns to p i the current value of clock. Event e i is the user event associated to point i: it represent a user interaction. Set P r represents the predecessors of i. Process U ser persistently chooses between launching or not an interactive event. CLOCK is a process that ticks forever. Process Score is parametrized by the graph in Fig. 1.2. We post a constraint 0 ≤ p i < n ∞ to allow that clock + 1 < p i be eventually deduced even in absence of user interactions.

iP
B,SL A,ER SC EC (0 , ∞ ) (0 , ∞ ) (0 , ∞ ) ∆ R ∆ G ∆ S ∆ C ∆ L SG,
CLOCK(k) def = tell (clock = k) next CLOCK(k + 1) Score g def = P oints lV,A Intervals lA U ser lV CLOCK(0) ! i∈V 0 ≤ p i < n ∞ 1.1.2. Discussion
We presented a model for interactive scores based upon ntcc to represent temporal reductions and interactive events. The model is close in spirit to the petri nets model. An advantage of petri nets is that transformation from interactive scores to petri nets is simpler, and synchronization is easier when it depends, for instance, on the transitions that precede a place. An advantage of ntcc is that it can easily represent global constraints such number of simultaneous temporal constraints and also temporal reductions.

The ntcc model could be implemented. We plan to implement our system using Ntccrt [TOR 09], a real-time capable interpreter for ntcc. There is an issue with the correctness of the implementation: it heavily relies on propagation. Each time an interactive point is launched, we add a constraint, and we propagate. Does propagation preserves the coherence of the model? Do we need to perform search or propagation is enough? This is an open issue. There is another open issue: We believe that the denotation of a score will help us to understand the behavior of the score without analyzing its operational semantics, to specify properties, and to prove its correctness.

Timed Conditional Branching Model

There is neither a formal model nor a special-purpose application to support conditional branching in interactive multimedia. Using conditional branching, a designer can create scenarios with choices. The user and the system can take decisions during performance with the degree of freedom described by the designer, and define when a loop ends; for instance, when the user changes the value of a variable.

In the domain of composition of interactive music, there are applications such as Ableton Live1 . Using Live, a composer can write loops and a musician can control different parameters of the piece during performance. Unfortunately, the means of interaction and the synchronization patterns are limited.

To express complex synchronization patterns and conditions, it was shown in [TOR 10] that interactive scores can describe temporal relations together with conditional branching. In this section, we recall the conditional-branching timed model of interactive scores based upon the Non-deterministic Timed Concurrent Constraint (ntcc) calculus. We explain how to model temporal relations, conditional branching and discrete interactive events in a single model. We also present performance results of a prototype and we discuss the advantages and limitations of the model.

Specification of the model

A score is defined by a tuple composed by a set of points and a set of intervals. A temporal object is a type of interval.

Points

We say that a point p is a predecessor of q if there is a relation p before q. Analogically, a point p is a successor of r if there is a relation r before p. A Point is defined by p = b p , b s , where b p and b s represent the behavior of the point. Behavior b p defines whether the point Waits until All its predecessors transfer the control to it or it only Wait for the First of them. Behavior b s defines whether the point CHooses one successor which condition holds to transfer the control to it, or it does Not CHoose, transferring the control to all its successors.

Intervals

An interval p before q means that the system waits a certain time to transfer the control (jumps) from p to q if the condition in the interval holds, and it executes a process throughout its duration. An interval is composed by a start point, an end point, a condition, a duration, an interpretation for the condition, a local constraint, a process, parameters for the process, children, and local variables.

There are two types of intervals. Timed conditional relations have a condition and an interpretation, but they do not have children, their local constraint is true, and their process is silence2 . Temporal objects may have children, local variables and a local constraint, but their condition is true, and their interpretation is "when the condition is true, transfer the control from the start point to the end point".

A timed conditional relation includes the start and end points involved in the relation and the condition that determines whether the control jumps from start to the end point. There are two possible values for the interpretation of the condition: (i) when means that if condition holds, the control jumps to the end point; and (ii) unless means that if the condition does not hold or its value cannot be deduced from the environment, the control jumps to end point. A temporal object is an interval where the start point launches a new instance of a temporal object and the end point finishes such instance. Its local variables and local constraint can be used by its children and process to synchronize each other. The points have the following behavior. The end point of C is enabled for choice, and the other points transfer the control to all their successors. The start point of C waits for all its predecessors and all the other points only wait for the first predecessor that transfers the control to them. Object A 's process is silence, its children are B, C and D and its local variable is f inish. Note that the silence between D and C lasts longer during execution because of the behavior of the points.

Ntcc model

We model points and intervals as processes in ntcc. The definition of an interval can be used for both timed conditional relations and temporal objects. We represent the score as a graph; however, there is not yet a procedure to remove zero-delay arcs. Predecessors and successors of point i are P redec i and Succ i , respectively. For simplicity, we do not include hierarchy, we only model the interpretation when, we can only declare a single interval between two points, and we can only execute a single instance of an interval at the same time.

Points

We only model points that choose among their successors (ChoiceP oint), points that transfer the control to all their successors (JumpT oAllP oint), and points that wait for all their predecessors to transfer the control to them (W aitF orAllP oint).

To know if at least one point has transferred the control to the point i, we ask to the store if j∈P oints Arrived(i, j) can be deduced from the store. When all the expected predecessors transfer the control to the point i, we add a constraint ActiveP oints i . Analogaly, when a point i transfers the control to a point j, we add the constraint ControlT ranf erred(j, i). In order to represent choice in the example of The following definition uses the parallel composition agent to transfer the control to all the successors of the point i.

T oAll i def = j∈P when Succs(i, j) do tell (ControlT ransf erred(j, i)))

tell (ActiveP oints i )

Using the definition T oAll i , we define the two types of point that transfer the control to all its successors. To wait for all the predecessors, we ask the store if the constraint Arrived = P redec holds.

JumpT oAllP oint i def = ! when j∈P oints Arrived(i, j) do T oAll i 1.2.2.2. Intervals Process I waits until at least one point transfers the control to its start point i, and at least one point has chosen to transfer the control to its destination j. Afterwards, it waits until the duration of the interval is over3 . Finally, it transfers the control from i to j. I i,j,d def = !(tell (P redec(j, i)) tell (Succ(i, j)))

! when k∈P ControlT ransf erred(j, k) ∧ k∈P Arrived(i, k) do( next d (tell(Arrived(j, i)) P redecessorsW ait(i, j)))

P redecessorsW ait adds the constraint Arrived(j, i) until the point j is active.

Example

We can represent the score in Fig. 1.3 in ntcc. The start point of textscc waits for all the point, the end of textscc chooses a point and the other points jump to all points.

The intervals have the duration described in Fig. 1.3. We also need to model a user making choices. U ser tells to the store that f inish is not true during the first n time units, then it tells that f inish is true. It is initialized with i = 0. An advantage of ntcc is that the constraint i ≥ n can be easily replaced by more complex ones; for instance, "there are only three active points at this moment in the score".

U ser n (i) def = when i ≥ n do tell (f inish) unless i ≥ n next tell (¬f inish) next U ser n (i + 1)

Results and Discussion

Performance results are described in [TOR 10]. They ran a prototype of the model over Ntccrt. The tests were performed on an iMac 2.6 GHz with 2 GB of RAM under Mac OS 10.5.7. They compiled it with GCC 4.2 and Gecode 3.2.2. The authors of the Continuator [PAC 02] argue that a multimedia interaction system with a response time less than 30 ms is able to interact in real-time with even a very fast guitar player. Response time was less than 30ms for the conditional branching model for up to 500 temporal objects.

An advantage of ntcc with respect to previous models of interactive scores, Pure Data (PD), Max and Petri Nets is representing declaratively conditions by the means of constraints. Complex conditions, in particular those with an unknown number of parameters, are difficult to model in Max or PD [PUC 98] . To model generic conditions in Max or PD, we would have to define each condition either in a new patch or in a predefined library. In Petri nets, we would have to define a net for each condition. A disadvantage of this model is that we cannot always synchronize two objects to happen at the same time: one reason is choice and the other is that using jumps is not always possible to respect the constraints on their durations.

Figure 1 . 1 .

 11 Figure 1.2 represents the graph associated to the score in Figure 1.1. It does not contain zero-delays.

Figure 1 . 2 .

 12 Figure 1.2. Graph representing the score in Fig. 1.1. Label St denotes the start point of t and Et the end point of t. Double arrows represent the duration of temporal objects. Simple arrows with no label represent the constraints imposed by the hierarchy, its duration is [0, ∞).

  Figure 1.3 describes a score with a loop. During the execution, the system plays a silence of one second. After the silence, it plays the sound B during three seconds and simultaneously it turns on the lights D for one second. After the sound B, it plays a silence of one second, then it plays video B. If the variable f inish becomes true, it ends the scenario after playing the video C; otherwise, it jumps back to the beginning of the first silence after playing the video C.

Figure 1

 1 Figure 1.3. A score with a user-controlled loop.

  Fig. 1.1, we use the constraint f inish. ChoiceP oint i,a,b def = ! when j∈P Arrived(i, j) do (tell (ActiveP oints i ) when f inish do tell (ControlT ransf erred(a, i)) +when ¬f inish do tell (ControlT ransf erred(b, i)))

http://www.ableton.com/live/

silence is a process that does nothing.

next d is a process next nested d times (next(next(next...).


Jean-Christophe C Audru

Bruno Capdeville

Jean-Jacques J Salindre

N Nérée

Éliane Mouly

Actions towards seismic risk

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The French Antilles National Earthquake Plan

>

  Initiated in the French Antilles in 2007 > To reduce the vulnerability of people and buildings > Budget ~100 M€ for Martinique actions > Actions dedicated to • public information and preparedness • professionals training courses • scientific knowledge • reinforcement and retrofit > The Réplik group of part ers = pu li ser i es, academy, regional and general council, foundations, experts, architects, mayors association, firemen, army, communication society, etc. > Annual information and preparedness campaign (all year + November) > Targets residents (adults, pupils) and tourists • Brochures, magnets, mouse pads, calendars, whistles, roadmaps, etc. • Exhibition stand at the annual Housing Fair (October) • Dedicated web site with a download area > Painting and sculpture exhibition, theatre skits, newspaper ads > TV, cinema and radio spots > Expanded partnership : pharmacies, medical offices, insurance firms, post offices, petrol outlets, commercial centers, TV channels etc. Specific information to employees and to construction professionals > Conferences for employees during working hours, Do and do 't flyers adapted to the industrial purposes (6000 people since 2008) > Training courses for construction professionals, resistant, possibly with some minor work to be performed • G2 : The retrofit/strengthening is fiscally sustainable • G3 : Reconstruction fiscally better than retrofit/strengthening > Current diagnosis of private schools Reinforcement & retrofit solutions > adapted to local schools, to local technical possibilities Reinforcement & retrofit projects > To day, reinforcement and reconstruction of 31 schools for 71 M€ > These are all subsidized operations (~ 80%): FPRNM (Special Natural Hazard Fund), FEDER (Europe), Government, Regional and General councils, ... my house in the state of the art > More training sessions for const. professionals > More retrofit projects > Diagnosis of governement buildings > New seismic codes applicable from next May > Updating of the Réplik internet site > New Réplik information

Tha k you for your atte tio a d … pa oli !


J Malick

F Roupin

On the

Keywords: Lagrangian duality, combinatorial optimization, 0-1 quadratic programming, nonlinear programming, semidefinite programming, branch-and-bound, quasi-Newton

bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds for 0-1 quadratic problems leading to quasi-Newton methods Jérôme Malick • Frédéric Roupin the date of receipt and acceptance should be inserted later This paper is dedicated to Claude Lemaréchal on the occasion of his 65th birthday. We take this opportunity to thank him deeply for the great moments we have had discussing with him (not only about math). His vision and his ability to put ideas into words has helped us deepen our understanding of optimization

study new bounds for combinatorial optimization, obtained by Lagrangian duality and combined with nonlinear optimization. To some extent, these bounds generalize and enhance the usual semidefinite programming bounds.

Semidefinite programming (SDP) is a prolific branch of convex optimization (see handbooks [START_REF] Saigal | Handbook of Semidefinite Programming[END_REF] and [START_REF] Anjos | Handbook of semidefinite, conic and polynomial optimization[END_REF]). Since the '90s, and in particular since the celebrated paper [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF] about max-cut, semidefinite programming has become a fundamental bounding technique for combinatorial optimization. The recipe to formulate SDP relaxations of combinatorial problems [START_REF] Poljak | A recipe for semidefinite relaxation for (0,1)-quadratic programming[END_REF] has an illuminating interpretation by Lagrangian duality, as explained by Claude Lemaréchal in the pedagogic paper [START_REF] Lemaréchal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF]. Recent research has spawned a focus on the power and the limitations of semidefinite programming for solving combinatorial optimization problems. The current sentiment is, roughly speaking, that semidefinite programming bounds have good tightness, but demand much time to be computed, which often prevents their direct use in branch-and-bound algorithms or other exact resolution schemes. To answer this challenge, a strategy is to use the problem structure, for example by exploiting sparsity (see e.g. [START_REF] Kim | Exploiting sparsity in sdp relaxation for sensor network localization[END_REF]), by reducing the size of problems by using hidden symmetries (see e.g. [START_REF] De Klerk | Exploiting group symmetry in semidefinite programming relaxations of the quadratic assignment problem[END_REF]), by developing new algorithmic schemes (see e.g. [START_REF] Malick | Regularization methods for semidefinite programming[END_REF]), or by adapting algorithms to special problem classes (see e.g. [START_REF] Rendl | Solving max-cut to optimality by intersecting semidefinite and polyedral relaxations[END_REF]).

In this paper, we use a different approach: we do not consider the standard SDP bounds; rather we construct new SDP-based bounds that are less tight (they deteriorate the usual SDP bounds) but are easier to compute (they are amenable to efficient standard nonlinear optimization techniques). The influence of Claude Lemaréchal is visible from the basic approach (Lagrangian duality [START_REF]Lagrangian relaxation[END_REF]) to the practical resolution (limited memory quasi-Newton method [START_REF] Gilbert | Some numerical experiments with variable-storage quasi-Newton algorithms[END_REF]); we also extensively use convex analysis [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms[END_REF].

1.2 A family of SDP bounds for 0-1 quadratic problems leading to a quasi-Newton method This article presents a family of semidefinite bounds for binary quadratic problems subject to linear or quadratic constraints, that is, a general problem of the form    min x ⊤ M 0 x + m 0 ⊤ x x ⊤ M j x + m j ⊤ x + c j = (or ≤) 0 j = 1, . . . , m x ∈ {0, 1} n

(1) with m j ∈ R n , c j ∈ R and symmetric matrices M j . This general form contains for example 0-1 quadratic problems with linear constraints (when M j = 0), and 0-1 linear problems (when moreover M 0 = 0), that are in particular models for various problems in operations research. Notice that the equality and inequality constraints in (1) can be both constraints arising when modeling the problem, or redundant constraints added to strengthen the bounds. To illustrate our approach, we consider in this paper three combinatorial optimization problems of the form (1): unconstrained 0-1 quadratic programming problems, heaviest subgraph problems, and graph bisection problems (see details and references in Section 2.2).

The objective of this article is to develop bounding techniques for (1), trading computing time for a deterioration of the quality of the usual semidefinite bounds, in order to enhance the efficiency of exact resolution schemes. We introduce a new family of SDP-based bounds constructed by duality after handily rewriting the rank-one constraint appearing in the SDP reformulation of the problems (see Section 2.4). The computation of these bounds relies on robust tools of linear algebra and nonlinear optimization (more precisely, computing the bounds amounts to an eigendecomposition of a matrix, and improving the tightness can be done with Newton-like algorithms; see Section 4.3). The bounds turn out to have good ratios of tightness to computing time (see Section 5), which is a key property in view of being used as a bounding procedure within branch-and-bound methods to solve 0-1 quadratic problems to optimality. Another important property is that we have a real parameter (α in the sequel) that controls the tightness level and the computing cost of these bounds. Some of the underlying ideas in this work are not completely new: they are present in [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF], and used in the recent paper [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF] to solve to optimality a special hard combinatorial optimization problem. The main goals of the present paper are the following:

to clarify and simplify the approach by introducing a family of semidefinite bounds with direct interpretation from Lagrangian duality; to give a general methodology to construct and to use these semidefinite bounds for general 0-1 quadratic programming; to make a numerical study of these bounds: illustrating their features, comparing them with the usual bounds (in terms of tightness, computing time, and the balance of both), and advocating that they have interesting qualities in view of being embedded within branch-and-bound algorithms.

Here is the outline of the paper. In Section 2, we introduce the notation, three problems of the form (1), and transformations to come up with an equivalent SDP formulation of the problem (see forthcoming equation ( 15)). In Section 3, we apply the mechanism of Lagrangian duality to the reformulation, which gives us the new bounds. We study in Section 4 the properties of this family of bounds and discuss how to use them in practice. Section 5 reports numerical tests of these bounds on various standard test-problems.

2 Semidefinite reformulations of 0-1 quadratic problems 2.1 Notation, projections in symmetric matrix spaces Let S p be the space of p × p symmetric matrices, equipped with the norm • associated with its usual inner product A, B := ∑ i, j A i j B i j = trace(AB). This inner product is relevant for our purposes through the relation x ⊤ A x = A, xx ⊤ for all x ∈ R p and A ∈ S p . In other words, one can transform quadratic functions on R p to linear functions on S p with respect to the rank-one symmetric matrix X = xx ⊤ . This will be used to obtain the forthcoming SDP formulation (15). Note also that the binary constraint x i ∈ {-1, 1} reads X ii = 1 on X = xx ⊤ . Introducing the operator diag : S p → R p and the vector of all ones e ∈ R p , this gives us the constraint diag(X) = e that appears in (15). For more details on SDP formulations of combinatorial problems, see [START_REF] Poljak | A recipe for semidefinite relaxation for (0,1)-quadratic programming[END_REF] and [START_REF] Lemaréchal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF].

The projections (with respect to • ) of a matrix A ∈ S n onto the closed convex cone S + p and its polar cone S - p are denoted, respectively, by

A + := argmin X 0 X -A and A -:= argmin X 0 X -A .
Moreau's theorem (see [HUL93, Chap. III]) states that the projections are characterized by

A + , A -= 0 and A = A + + A -. (2)
In fact, it is well-known (see e.g. [START_REF] Higham | Computing a nearest symmetric positive semidefinite matrix[END_REF]) that we have an explicit formula for these projections as follows. Let A = ∑ i λ i p i p i ⊤ be a spectral decomposition of A with eigenvalues λ i and eigenvectors p i , assumed to be pairwise orthogonal and of length one; then

A + = ∑ λ i >0 λ i p i p i ⊤ and A -= ∑ λ i <0 λ i p i p i ⊤ .
Note finally that we have for A ∈ S p (-A) -= -A + .

(3)

Examples of binary quadratic optimization problems

The developments of this paper apply to general binary quadratic problems of the form (1). We illustrate them on the following three classical combinatorial optimization problems.

Example 1 (Unconstrained quadratic 0-1 problem) Consider the problem of minimizing a quadratic function with binary variables

min x ⊤ M 0 x + m 0 ⊤ x x ∈ {0, 1} n . ( 4
)
This problem is known to be NP-hard in general, and to have numerous applications (see e.g. an application in medical field in [PIS + 04]). We also mention the celebrated max-cut problem on a (undirected) graph (see e.g. [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF]) which can be formulated as

min x ⊤ L x x ∈ {-1, 1} n (5)
with L the Laplacian matrix of the graph. It is well-known that these two problems (4) and (5) are equivalent. In particular, (4) can be reformulated as a max-cut on a graph with unrestricted weights.

⊓ ⊔

Example 2 (Heaviest subgraph) Consider an undirected weighted graph G = (V, E) with n vertices, and denote by W = (w i j ) i j its weight-matrix (which is a symmetric n × n matrix).

For an integer k in {1, . . . , n}, the heaviest k-subgraph problem (also called k-cluster problem; see references in [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF]) consists in determining a subset S of k vertices such that the total edge weight of the subgraph induced by S is maximized. This problem is directly modeled in the form (1) max

x ⊤ W x e ⊤ x = k, x ∈ {0, 1} n . (6)
In view of having better bounds when relaxing, we introduce the n so-called "product constraints" (standard reinforcement constraints, see e.g. [START_REF] Sherali | A hierarchy of relaxations between the continuous and convex hull representations for zero-one programming problems[END_REF], [START_REF] Lovász | Cones of matrices and set-functions and 0-1 optimization[END_REF])

(e ⊤ x -k) x i = 0, for i = 1, . . . , n. (7)
In some sense, adding these n equality constraints is optimal (see e.g. [START_REF] Lemaréchal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF], [START_REF] Faye | Partial lagrangian for general quadratic programming. 4'OR[END_REF]). This leads again to a problem of the form (1)

   max x ⊤ W x e ⊤ x = k, x ∈ {0, 1} n , x ⊤ A i x -k x i = 0, for i = 1, . . . , n (8)
with the appropriate matrices A i ∈ S n (see [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF]Eq. (3)] for more details).

⊓ ⊔

Example 3 (Graph bisection) With the notation of the previous example, the graph bisection problem consists in dividing the nodes of G into two subsets of k and nk vertices such that the total weights of edges that have end-points in different sets is minimal. This problem can be formulated as min

x ⊤ W x e ⊤ x = 2k -n, x ∈ {-1, 1} n . (9)
Similarly to heaviest subgraph problem (8), the constraints can be strengthened as

   min x ⊤ W x e ⊤ x = 2k -n, x ∈ {-1, 1} n , x ⊤ A i x -(2k -n) x i = 0, for i = 1, . . . , n. ( 10
)
Moreover this problem can be obviously written as (1) after a change of variables. SDP techniques have been used for solving this problem to optimality (see [START_REF] Karisch | Solving graph bisection problems with semidefinite programming[END_REF]). ⊓ ⊔ 2.3 Reformulation with a rank-one constraint (semidefinite lifting)

We apply the standard lifting up to the symmetric matrix space by equivalently writing the 0-1 quadratic problem (1) as a linear SDP problem with a rank-one constraint (see more in [START_REF] Saigal | Handbook of Semidefinite Programming[END_REF]). Though all the development can be made with {0, 1}-variables, it is simpler to work with {-1, 1} variables when introducing the "spherical constraint" (see next section).

There are two equivalent ways to proceed when changing variables and lifting to SDP:

1. Lift and change. First, apply the "recipe" (see e.g. [START_REF] Poljak | A recipe for semidefinite relaxation for (0,1)-quadratic programming[END_REF] or [START_REF] Lemaréchal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF]) to transform the 0-1 quadratic problem (1) to a SDP problem with respect to

Y = 1 x ⊤ x X ∈ S n+1 .
Second, apply the change of variables X ← UYU ⊤ given by the block-matrices

U := 1 0 1 2 e 1 2 I n and U -1 = 1 0 -e 2I n .
2. Change and lift. We can also do it the other way around. First, apply (if necessary) the change of variable {0, 1} → {-1, 1} to (1). In fact, the transformation [1; x] = U[1; y] with the above U can be used to make the change of variables and to homogenize the quadratic forms q j (x) = x ⊤ A j x + b j ⊤ x + c j at the same time. Second, lift the problem to the matrix space with a semidefinite rank-one matrix of the form X = xx ⊤ ∈ S n+1 .

The two above transformations are known to be equivalent (see e.g. [START_REF] Helmberg | Combining semidefinite and polyhedral relaxations for integer programs[END_REF], [START_REF] Helmberg | Semidefinite Programming for Combinatorial Optimization[END_REF]). Note that there are several degrees of freedom in the reformulation (although all the operations can be done automatically). For example, the constants c j can be treated in different ways. The formulation choices may have an impact on the behavior of solvers.

Finally we get the following linear SDP problem with respect to

X ∈ S n+1    min Q 0 , X Q j , X = (or ≤) q j j = 1, . . . , m diag(X) = e, X 0, rank(X) = 1 (11)
with γ j , q j ∈ R such that γ jq j = c j , and with the (n + 1) × (n + 1)-matrices Q j defined by

Q j := U ⊤ γ j 1 2 m j ⊤ 1 2 m j M j U.
For example, problem (4) has the equivalent SDP formulation

   min Q, X diag(X) = e X 0, rank(X) = 1 with Q := 1 4 e ⊤ M 0 e + 2m 0 ⊤ e e ⊤ M 0 + m 0 ⊤ M 0 e + m 0 M 0 . ( 12
)
The combinatorial difficulty of the initial problem (1) is now concentrated in the rankone constraint of (11). The SDP relaxation of this problem consists of dropping the rank-one constraint; the resulting problem is thus the linear SDP problem

   min Q 0 , X Q j , X = (or ≤) q j j = 1, . . . , m diag(X) = e, X 0. ( 13
)
In this paper, we do not drop the rank constraint, and we follow instead the idea of [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF] to write this constraint in a more handy way, as a norm constraint.

2.4 Reformulation of the rank-one constraint (as a spherical constraint)

An interesting property [Mal07, Th. 1] is that, for X 0 with diag(X) = e, we have

X ≤ n + 1 and moreover X = n + 1 ⇐⇒ rank X = 1 . ( 14
)
Therefore we can replace the rank-one constraint in formulation (11) by the equivalent constraint X 2 = (n + 1) 2 , called the "spherical constraint" in [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF]. Our problem (1) can thus be equivalently written as:

   min Q 0 , X Q j , X = (or ≤) q j j = 1, . . . , m diag(X) = e, X 0, X 2 = (n + 1) 2 . (15)
In other words, we have transformed a quadratic problem in R n+1 with quadratic constraints into a linear SDP problem with one single quadratic constraint. This quadratic constraint carries the nonconvexity and the combinatorial difficulty of the initial problem. Obviously, dropping it take us back to the SDP relaxation (13).

Semidefinite bounds by Lagrangian duality

In this section, we attack the last reformulation (15) of the quadratic binary problem (1), by Lagrangian duality. We will apply, simply and directly, the usual duality mechanism (see e.g. [BV04, Chap. 5] and [HUL93, Chap. XII]). To do so, we separate the m affine constraints of (15) into m I inequalities and m E equalities (with m I + m E = m, and m E ≥ n since there are already the n constraints X ii = 1). So we work from now on with

       min Q, X A i , X ≤ a i , i = 1, . . . , m I B i , X = b i , i = 1, . . . , m E X 0, X 2 = (n + 1) 2 (16)
which is (15) with new notation. We emphasize that there has not been any relaxation so far; in other words, the optimal value of ( 16) is equal to the one of the original problem (1): val (16) = val (1).

We introduce the dual space

D := R m I + × R m E × S - n+1 × R
, and define the Lagrangian, a function of the primal variable X ∈ S n+1 and the dual variables

(λ , µ, Z, α) ∈ D L(X; λ , µ, Z, α) = Q, X + m I ∑ i=1 λ i ( A i , X -a i ) + m E ∑ i=1 µ i ( B i , X -b i ) (17) + α 2 ( X 2 -(n + 1) 2 ) + Z, X .
The associated dual function is the concave function defined as

θ (λ , µ, Z, α) := inf X∈S n+1 L(X; λ , µ, Z, α). ( 18
)
We already know by the weak duality result ([HUL93, XII.2.1.5]), that each value of θ gives a lower bound on the optimal value of (16); in other words, we have for all

(λ , µ, Z, α) ∈ D θ (λ , µ, Z, α) ≤ val (16). ( 19
)
The best of these bounds is the optimal value of the associated dual problem sup

(λ ,µ,Z,α)∈D θ (λ , µ, Z, α). (20)
So far, we have just given definitions and general properties; the approach has a real interest only if we can compute (some of) these bounds. In general, solving the optimization problem (18) -and therefore (20) a fortiori -is not easy as it has a nonlinear objective function and a semidefinite conic constraint. In fact, the computability of this problem essentially relies on the sign of α, as follows.

First of all, when α < 0, observe that θ (λ , µ, Z, α) = -∞ so that this "bound" is useless. Second, the case when α = 0 deserves a special treatment. Indeed note that the Lagrangian L(X; λ , µ, Z, 0) corresponds to the standard Lagrangian for linear SDP programming. Thus we have an explicit expression of θ (λ , µ, Z, 0) that gives us the usual linear SDP dual problem (which corresponds to the dual of (13))

   max -a ⊤ λ -b ⊤ µ A(λ ) + B(µ) + Q + Z = 0 Z 0 (21)
where we set

A(λ ) = m I ∑ i=1 λ i A i and B(µ) = m E ∑ i=1 µ i B i .
In other words, the best of the bounds θ (λ , µ, Z, 0) for the dual variables (λ , µ, Z) is simply the usual SDP bound of (1). There exist many efficient algorithms and software to compute this SDP bound (see e.g. [START_REF] Helmberg | An interior point method for semidefinite programming[END_REF], [HR00]). In our numerical experiments, we use the two software packages [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF] and [START_REF] Helmberg | A C++ implementation of the Spectral Bundle Method[END_REF] known to be among the best ones. Finally, problem (18) with α > 0 can be solved explicitly, as shown by the next two results. We denote the positive real numbers by R ++ .

Theorem 1 (Dual function) Given dual variables (λ , µ, Z, α) ∈ D with α > 0, the dual function can be written

θ (λ , µ, Z, α) = -a ⊤ λ -b ⊤ µ - 1 2α Q + A(λ ) + B(µ) + Z 2 - α 2 (n + 1) 2 .
This expression is differentiable with respect to

(λ , µ, Z, α) ∈ R m I × R m E × S n+1 × R ++ ,
and we know the partial differentials; in particular:

∂ λ i θ (λ , µ, Z, α) = - 1 α A i , Q + A(λ ) + B(µ) + Z -a i ∂ µ i θ (λ , µ, Z, α) = - 1 α B i , Q + A(λ ) + B(µ) + Z -b i .
Proof Note that the Lagrangian (17) can be written as

L(X; λ , µ, Z, α) = α X 2 /2 + X, Q 1 (λ , µ, Z) + q(λ , µ, α)
with the help of

Q 1 (λ , µ, Z) := Q + A(λ ) + B(µ) + Z, and q(λ , µ, α) := - α 2 (n + 1) 2 -a ⊤ λ -b ⊤ µ.
As a function of X, the Lagrangian is obviously strongly convex and differentiable; thus it admits a unique minimizer characterized by 0

= ∇ X L(X; λ , µ, Z, α) = αX + Q 1 (λ , µ, Z). This gives θ (λ , µ, Z, α) = - 1 2α Q 1 (λ , µ, Z) 2 + q(λ , µ, α) (22)
which is the form we wanted. This function is clearly differentiable with respect to all the dual variables; the expressions of the differentials follow from standard calculus rules. ⊓ ⊔ Theorem 2 (Simplified dual function) Given dual variables (λ , µ, α) with α > 0, the dual function can be maximized over Z: the simplified dual function is

Θ (λ , µ, α) := max Z 0 θ (λ , µ, Z, α) = -a ⊤ λ -b ⊤ µ - 1 2α (Q + A(λ ) + B(µ)) - 2 - α 2 (n + 1) 2 .

This expression is differentiable at any

(λ , µ, α) ∈ R m I × R m E × R ++ ,

and we have

∂ λ i Θ (λ , µ, α) = - 1 α A i , (Q + A(λ ) + B(µ)) --a i ∂ µ i Θ (λ , µ, α) = - 1 α B i , (Q + A(λ ) + B(µ)) --b i . Proof Isolate matrix Q 2 (λ , µ) = Q + A(λ ) + B(µ) in the expression of θ (λ , µ, Z, α) of Theorem 1. Observe that the solution of min Z 0 Q 2 (λ , µ) + Z is the projection (-Q 2 (λ , µ)) -= -Q 2 (λ , µ) + by (3). Therefore the minimum value is then -Q 2 (λ , µ) + + Q 2 (λ , µ) = Q 2 (λ , µ) -by (2).
With the expression of θ (λ , µ, Z, α) of Theorem 1, this gives the expression of Θ (λ , µ, α).

To differentiate Θ , we want to use the theorem of differentiability of a max (or min) function (see e.g. [HUL93, VI.4.4.5]) which asks for compactness of the "index set" on which the maximum is taken. Fix any ( λ , μ, ᾱ) ∈ R m I × R m E × R ++ and consider a closed ball V around them. Observe now that the set U = -(Q 2 (V )) + is compact as well (since Q 2 and the projection are both continuous) and that we can write

Θ (λ , µ, α) = max Z∈U∩S - n+1 θ (λ , µ, Z, α) for all (λ , µ, α) ∈ V.
Since the maximum is now taken on a compact set, [HUL93, VI.4.4.5] gives the differentiability around ( λ , μ, ᾱ) together with the expression of partial differentials.

Let us make the calculation of the partial differential with respect to λ i explicit. By [HUL93, VI.4.4.5], we have

∂ λ i Θ (λ , µ, α) = ∂ λ i θ (λ , µ, -Q 2 (λ , µ) + , α)
. The expression of ∂ λ i θ of Theorem 1 yields, with the help of (2)

∂ λ i Θ (λ , µ, α) = - 1 α A i , Q 2 (λ , µ) -Q 2 (λ , µ) + -a i = 1 α A i , Q 2 (λ , µ) --a i
We obtain the other partial differentials the same way.

⊓ ⊔

This theorem shows that for any dual variables (λ , µ, α) with α > 0 we have an explicit expression of the dual lower bound Θ (λ , µ, α). In the remainder of this paper, we study these new bounds -theoretically and numerically.

Study of the family of semidefinite bounds

Consider the bounds for the 0-1 quadratic problem (1) formulated as ( 16), introduced by duality in the previous section; we study in this section the family of bounds

Θ (λ , µ, α) = -a ⊤ λ -b ⊤ µ - 1 2α (Q + A(λ ) + B(µ)) - 2 - α 2 (n + 1) 2 (23)
with respect to (λ , µ, α) ∈ D + where

D + := R m I + × R m E × R ++ .
Obviously, these are indeed lower bounds of the optimal value of (1), since in view of (19) and Theorem 2 we have

Θ (λ , µ, α) ≤ val (16) = val (1). ( 24
)
We study in this section some important properties of the bound Θ (λ , µ, α) (computation cost in Section 4.1, tightness in Section 4.2) and we discuss how to use them in practice within a branch-and-bound procedure (Section 4.3). A numerical study on our three testproblems is done afterwards in Section 5.

Computing of the bounds

The bounds Θ (λ , µ, α) are computed by standard linear algebra:

1. The initial quadratic problem is reformulated as ( 16) by a sequence of elementary operations that could be fully automatized (to be implemented in SDP S [START_REF] Roupin | From linear to semidefinite programming: an algorithm to obtain semidefinite relaxations for bivalent quadratic problems[END_REF]). 2. The main operation is then to assemble the matrix Q + A(λ ) + B(µ) and to project it onto the cone of semidefinite negative matrices S - n+1 .

As recalled in Section 2.1, projecting onto S - n+1 is done by computing an eigendecomposition of the matrix to project. The dominant cost of the computation of Θ (λ , µ, α) is thus of order O(n 3 ). This may prevent the use of these bounds for large graphs. Since we aim at developing a bouding technique to be embedded within a branch-and-bound algorithm, this is not a limitation: regarding exact resolution of our three test-problems, medium-size problems (n ≤ 500) are already challenging instances, to say the least. In this case, the eigendecomposition can be performed quickly by efficient linear algebra routines. Thus, the new bound are not too expensive in view of our targeted use; but, of course, they are of relevance only if they are not too weak. This is studied in the next section.

Before moving to tigthness, let us make the bounds for one of our test-problems explicit for illustration. This will also highlight the connection with [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF].

Example 4 (Bounds for the heaviest subgraph problem) Particularizing bound (23) for problem (8) yields, for µ

= (µ 1 , µ 2 ) ∈ R n+1 × R n+1 and α > 0, Θ heav (µ, α) := -b ⊤ µ - 1 2α (Q + A heav (µ)) - 2 - α 2 (n + 1) 2 ,
where we set

Q = U ⊤ 0 0 0 -W U and A heav (µ) := Diag(µ 1 ) +U ⊤    0 k 2 µ 2 ⊤ k 2 µ 2 - n+1 ∑ i=1 (µ 2 ) i A i   U.
We note that the above bounds have already appeared in [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF]: more specifically, these are not the bounds introduced in [MR11, Sec. 3], but they are essentially the same as the ones used in practice, implicitly given by [MR11, Lemma 3] as

Θ ′ heav (µ, α) := α -b ⊤ µ + 1 2 (-Q/α + A heav (µ)) + 2 - 1 2 (n + 1) 2 .
Indeed we have Θ heav (µ, α) = Θ ′ heav (-µ/α, α). In addition to the direct interpretation by duality, one could see a numerical interest in using the formulation Θ heav instead of Θ ′ heav to compute the bound. To get Θ heav indeed, we compute the eigendecomposition of a matrix whose entries are of the same scale as the data of the problem, whereas the entries of the matrix in Θ ′ heav are badly scaled when α gets small. ⊓ ⊔

How tight are the bounds?

This section shows that the new SDP bounds (23) are always less tight than the usual SDP bounds, but that one can approximate them in arbitrary precision. The first theorem goes one step further than (24) by showing that the new SDP bounds are in fact lower bounds on the usual SDP bound of the initial problem.

Theorem 3 (Less tight bounds) For any (λ , µ, α) ∈ D + , the bound Θ (λ , µ, α) is less tight than the usual SDP bound:

Θ (λ , µ, α) ≤ val (13) ≤ val (1).
Proof Denote by C the constraint set of the SDP relaxation (13), that then becomes

min Q, X X ∈ C . (25)
Recalling that ( 16) is just (15) written with other notation, observe that the feasible set of ( 16) is included in C , so that val (25) ≤ val (16) = val (1). ( 26)

Consider now α > 0 fixed, and gather together the terms Q, X and α( X 2 -(n + 1) 2 )/2 in (17) and (18): we then observe that

max θ (λ , µ, Z, α) (λ , µ, Z) ∈ R m I + × R m E × S - n+1 , (27)
is the dual problem of

min Q, X + α 2 ( X 2 -(n + 1) 2 ) X ∈ C . ( 28
)
By weak duality we have val (27) ≤ val (28). Since any X ∈ C is positive semidefinite with ones on the diagonal, (14) gives α( X 2 -(n + 1) 2 ) ≤ 0, so that we have furthermore val (28) ≤ val (25). Putting together the inequalities, we have

Θ (λ , µ, α) ≤ val (27) ≤ val (28) ≤ val (25) = val (13) ,
which allows us to conclude with (26).

⊓ ⊔

The next result somehow compensates the bad news of the previous one: it states that we can get close to the SDP bound with Θ (λ , µ, α). To prove the result, we need existence of solutions to (21), so we make a Slater-type assumption.

Theorem 4 (Arbitrarily tight bounds) Set Λ := R m I + × R m E , and assume that there exists a matrix X ∈ S n+1 such that

X ≻ 0, A i , X ≤ a i , i = 1, . . . , m I and B j , X = b j , j = 1, . . . , m E ( 29
)
Then the bounds Θ (λ , µ, α) satisfy the following two properties:

(i) for 0 < α ≤ β , max (λ ,µ)∈Λ Θ (λ , µ, β ) ≤ max (λ ,µ)∈Λ Θ (λ , µ, α) ≤ val (21) ;
(ii) the upper bound is reached when α > 0 vanishes:

lim α→0 max (λ ,µ)∈Λ Θ (λ , µ, α) = val (21).
Proof Using a standard result of convex optimization (apply e.g. [HUL93, VII.2.2.5] with c 0 (X) = λ max (-X)), the weak Slater assumption (29) gives the existence of dual solutions (i.e. solutions to (21) in our case), and together with existence of primal solution (the primal feasible set is compact), it also gives no duality gap (i.e. val (21) = val (13)). We use now the notation of the proof of Theorem 3 and we introduce the function

T (α) := max (λ ,µ)∈Λ Θ (λ , µ, α) = max (λ ,µ,Z)∈R m I + ×R m E ×S - n+1 θ (λ , µ, Z, α)
Note first that the Slater point of ( 13) is also a Slater point of (28) (that has the same constraint set, denoted C ). Therefore, the optimal values of ( 27) and (28) coincide, in other words, T (α) = val (28). Look now at (28): the constraint set C is compact (it is closed by definition and bounded by ( 14)), and the objective function is strongly convex (since α > 0). Therefore, there exists a unique solution X * ∈ C , and we can use the theorem of differentiability of a max function (see e.g. [HUL93, VI.4.4.5]): T is differentiable and

T ′ (α) = ( X * 2 -(n + 1) 2 )/2 ≤ 0.
The above negativity comes from (14) since matrices lying in C are in particular positive semidefinite with ones on the diagonal. Thus we have that T is nonincreasing on R ++ , and this proves (i) (since the last inequality in (i) comes from Theorem 3). We now turn to (ii). We denote by (λ * , µ * , Z * ) a solution to (21) (that exists by assumption); there holds θ (λ * , µ * , Z * , 0) = val (21). On the other hand, we have for α > 0

θ (λ * , µ * , Z * , α) ≤ Θ (λ * , µ * , α) ≤ max (λ ,µ)∈Λ Θ (λ , µ, α) ≤ val (21). ( 30
)
by definition of Θ and Theorem 3. Recall that θ is concave and upper-semicontinuous since it is defined in (18) as an infinum of linear functions L(X; •, •, •) (see e.g. [HUL93, IV.2.1.2]). Moreover, the domain of θ contains D + as well as the point (λ * , µ * , Z * , 0). We deduce that θ is continuous when α > 0 tends to zero; more specifically, we have

lim α→0 θ (λ * , µ * , Z * , α) = θ (λ * , µ * , Z * , 0) = val (21).
This limit together with (30) allows us to conclude the proof.

⊓ ⊔

Experiments in Section 5 illustrate this result: the optimal value of ( 21) is well approximated by taking α small (between 10 -2 and 10 -4 depending on the problem). Let us also mention that this convergence is accentuated in practice by rounding. The optimal value of 0-1 problem (1) is integer if the entries of the M j and m j are integers as well, so that rounding up Θ (λ , µ, α) still gives a lower bound to (1). Rounding up the SDP bound and Θ (λ , µ, α) for very small α would give the same bound (except in pathological cases).

Individually, the lower bound Θ (λ , µ, α) could have good tightness, but for a given triple (λ , µ, α), we still have no control on the quality of the bound. We explain in the next section how to use these bounds in practice.

How to use the bounds in practice?

Maximizing bounds. We note first that Theorem 4 yields another way to compute the SDP bound, specifically by solving max (λ ,µ,α)∈D + Θ (λ , µ, α) = val (21). However we do not use directly this property because the constraint α > 0 in D + complicates the picture: the constraint set is not closed, and handling it would require extra technicalities. We proceed otherwise, as follows.

Property (i) of Theorem 4 suggests that α acts like a control parameter to set the level of tightness. Fixing α > 0 yields a simpler optimization problem max

λ ∈R m I + , µ∈R m E Θ (λ , µ, α) (31)
with a differentiable objective function (by Theorem 2) and just nonnegativity constraints on the dual variable λ . This problem can be solved by any algorithm of nonlinear optimization (handling box constraints).

Newton-like algorithm. Though Θ is differentiable, it is not twice differentiable due to the projection onto the cone of negative semidefinite matrices in (23). Direct application of the Newton method is not possible, so we turn to Newton-like methods: quasi-Newton methods (see e.g. [START_REF] Bonnans | Numerical Optimization[END_REF]), or (inexact) semismooth Newton methods (see [START_REF] Qi | A nonsmooth version of Newton's method[END_REF] and [START_REF] Nocedal | Numerical Optimization[END_REF]).

In our numerical experiments, we use the limited memory quasi-Newton of [START_REF] Gilbert | Some numerical experiments with variable-storage quasi-Newton algorithms[END_REF] to solve (31). It proves to have good numerical behavior, as shown in the next section. Using a limited-memory version allows us to deal with problems (1) with possibly many constraints. In particular, one could add without concern redundant constraints to sharpen the bounds, as is done with the product-constraints in (8).

The parameter α controls the tightness level of the bounds. We also notice that (31) can be solved to various levels of accuracy, so that the stopping criterion of the quasi-Newton solver is a secondary practical parameter impacting the quality of the bound. In Section 5, the stopping criterion is set from 1e-4 to 1e-7 depending on the type of problem and on the value of α. We also mention that we scale the data initially to have better control on the stopping criterion: the constraints B i , X = b i are scaled so that B i are of norm 1.

We finally underline an important point. Throughout the run of the chosen maximizing algorithm, the current objective value Θ (λ k , µ k , α) is obviously a bound. So in practice within branch-and-bound, the run is stopped when the current value gets above the threshold given by the best known solution of (1). Only few iterations are often necessary to be able to prune the branch-and-bound tree (see the experiments with branch-and-bound in [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF] using Θ ′ heav of example 4).

Choosing the tightness parameter α. The question reduces now to how to choose the tightness level α. On the one hand, reducing α enhances the final gap obtained by the maximization (Theorem 4 (i)). On the other hand, reducing α also yields numerical difficulties. We can guess by the expression of Θ and its gradient (Theorem 2) that very small values of α can cause ill-conditioning. Numerical experiments confirm that solving the problem becomes more difficult as α decreases (see Section 5.2).

The goal now is to adapt α to deteriorate the tightness of the SDP bound in order to save computing time. Thus a balance has to be found when reducing α. In practice, one could do preliminary testing (as in Section 5.2) to choose an appropriate α for each problem.

Computational study and comparison

We illustrate the properties of the family of bounds Θ (λ , µ, α) on the three examples of Section 2.2. We denote the bounds (23) respectively by -Θ uqp for the unconstrained 0-1 quadratic programming problem (4), -Θ heav for the heaviest subgraph problem (2), -Θ bis for the graph bisection problem (3). Section 5.2 shows the profiles of the improvement in tightness when using α to control tightness. Afterwards, numerical experiments of Section 5.3 compare these new bounds with the standard semidefinite programming ones.

Machine, criterion and data

The numerical experiments have been carried out on a Pentium IV 2.2 GHz with 1 GB of RAM under Linux. The reported computing times are in seconds. In section 5.2, we evaluate the tightness of the bounds by the relative gap between the bound rounded up to the next integer, and the best known value of the initial problem. The best known solution can be either the optimal (computed by an exact resolution scheme), or suboptimal (computed by heuristics, when all exact schemes failed); more details are given for each case. Unconstrained 0-1 quadratic programming. We consider some of the biggest and most difficult problems of [START_REF] Billionnet | Using a mixed integer quadratic programming solver for the unconstrained quadratic 0-1 problem[END_REF]: randomly generated instances created with the generator of [START_REF] Pardalos | Computational aspects of a branch and bound algorithm for quadratic zero-one programming[END_REF] where n is the number of variables and d is the density of matrix M 0 in (4). For each couple of parameters among n = 100, 150, 200 and d = 30%, 80%, 100%

we have 10 instances of those problems. In tables 1 and 4, each reported result is thus the mean of the 10 instances for given n and d. The solution of (4) used to evaluate the gap for different bounds is as in [START_REF] Billionnet | Using a mixed integer quadratic programming solver for the unconstrained quadratic 0-1 problem[END_REF]: it is the exact solution for n = 100, and the best feasible available feasible solution for n = 150, 200.

Heaviest subgraph problems. We consider randomly generated instances used by several papers [START_REF] Billionnet | Different formulations for solving the heaviest k-subgraph problem[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] and [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF]. The parameters of these instances are: the size of the graph n, the graph density d, and the value of k. For each combination of parameters among: n = 80, 100, 300 d = 25%, 50%, 75% and k = n/4, n/2, 3n/4, we have 5 instances of problems. In tables 2 and 5, each reported result is the mean of the 15 instances for given n and d. For n = 80 and n = 100 the best known solutions are optimal (taken from [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF]). For n = 300, we used the heuristic of [MR11] (a greedy algorithm followed by a "two-opt" procedure) to generate a good feasible solution. Thus the real gaps between the bounds Θ heav and the optimal value are certainly smaller.

Graph bisection problem.

We consider some of the problems of [START_REF] Karisch | Solving graph bisection problems with semidefinite programming[END_REF]. We have 8 graphs of size n between 84 to 204; two of them are the weighted graphs ex84f and ex132f, the others (denoted ex*a) are unweighted. We get 32 instances of the bisection problem taking k = ⌊n/2⌋, ⌊3n/4⌋, ⌊7n/12⌋, ⌊13n/24 + 1/2⌋.

Note that for k = ⌊n/2⌋, we obtain the well-studied equicut problem. In the forthcoming tables 3 and 6, each reported result corresponds to one instance. The best known solutions are taken from [START_REF] Karisch | Solving graph bisection problems with semidefinite programming[END_REF], except for n ≥ 156. These solutions are optimal when n = 84 and are at most at 1% for n = 108 and n = 132. For n ≥ 156, we compute a suboptimal solution with a simple heuristic ("two-opt" procedure from randomly generated bisections).

Comparison of the bounds; role of α as a control parameter

We test and analyze the impact in terms of quality of the bound and computing time required for different choices of α. Such preliminary testing allows one to adjust α for each problem type. The three tables of this section report the time to optimize the bounds for a given value of α and the final gap between the bound and the best known solution of the initial combinatorial problem. Table 1 Numerical results for the unconstrained quadratic problem (4) (each entry is the mean of 10 results). The table reports the computing time of Θ uqp and gap(%) the averaged relative difference between the computed bounds and the best known solution of (4).

Θ uqp with α = 5.10 -2 Θ uqp α = 10 -2 n d time (s) gap(%) time (s) gap(%) time (s) gap(%
For unconstrained 0-1 quadratic programming, six different values of α are compared in table 1. Figure 1 gives an illustration of the convergence for a particular instance with n = 150 and d = 0.3. We note that using α = 10 -4 is more expensive than α = 10 -3 for only a small improvement of the gap.

For heaviest subgraph problems and for bisection problems, the picture is different. Tables 2 and 3 shows that the tightness greatly increases when α vanishes. In particular, taking α = 10 -4 provides a SDP-like bound in both cases. The computing cost differs though: for the heaviest subgraph problems, this is obtained within a reasonable computing time (when n ≤ 100), whereas it is more costly for graph bisection problems of the same sizes. More generally, we shall see in Section 5.3, that our best bounds are SDP-like for the three problems. Thus the corresponding gaps in Tables 1, 2, and 3 are almost equal to the SDP ones.

Problem

Θ heav with α = 10 -2 Θ heav with 5.10 -3 Θ heav with 10 Table 3 Numerical results for graph bisection problems. The table reports the computing time of Θ bis and gap(%), the averaged relative difference between the computed bounds and the best known solution of (8).

Comparison with standard semidefinite bounds

This section gives numerical comparisons of the new bounds and the SDP bound (21), in terms of tightness, computation time and the balance between both. To compute the SDP bound, we use two SDP solvers known to be efficient:

1. SB [START_REF] Helmberg | A C++ implementation of the Spectral Bundle Method[END_REF], a SDP solver based on the spectral bundle method [START_REF] Helmberg | A spectral bundle method for semidefinite programming[END_REF]. This software can handle large-scale problems and is known to be efficient in the context of combinatorial optimization. We use it with the default settings, except that: we set the stopping criterion to 1e-4 (instead of 1e-5), we activate (with -sh 1 and -si 1) the heuristic for finding a good starting point and the initial rescaling of the constraints. 2. CSDP [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF], a robust and efficient interior point solver. We use it with the default settings, except that we activate "Fastmode" and we set the stopping criterion objtol to 1e-5 (instead of 1e-8).

Remark 1 (Different formulations of the same SDP bound) Computing times depend obviously on the solvers, and in turn the performances of the solvers depend on the formulation of the bounds. There is no ambiguity in the SDP relaxation of (4), but there is some ambiguity in those of ( 8) and (10). It turns out that adding the n product constraints ( 7) is equivalent to adding the squared constraint (e ⊤ xk) 2 = 0 (see e.g. [START_REF] Lemaréchal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF]); this leads to two equivalent SDP relaxations. For the k-heaviest problem, preliminary tests of [START_REF] Malick | Numerical study of semidefinite bounds for the k-cluster problem[END_REF] showed that the formulation with the product constraints is best solved by SB, while the one with the squared constraint is best solved by CSDP. We have observed a similar behavior for the graph bisection problems. In the numerical experiments below, we take the most advantageous formulation for each solver. We emphasize this idea which is unusual (as far as we are aware) but essential to make fair comparisons. We see this idea as a secondary but important contribution of this paper. Tables 4, 5 and 6 report the computation times of the three solvers. Recall that the solvers do not compute the same bounds: to compare the SDP bound (21) and the maximum of Θ (µ, α) for fixed α, we report:

-gap% which is the relative difference between the two bounds, the computing times for SB and CSDP to converge, the computing times for SB and CSDP to achieve the final Θ . We emphasize that "gap" does not have the same meaning as in the previous tables: it is now the difference between Θ (µ k , α) and the SDP bound (21). We also note that to have a fair comparison, the bounds are not rounded up to the next integer as was done previously. Finally we mention that table 5 aggregates the results of [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF].

The three tables show that the new bounds provide SDP-quality bounds more quickly. More specifically, the difference between the final Θ and the SDP bound is always smaller than 0.1% (recall moreover that in practice it would be reduced again by rounding). Though SB and CSDP eventually converge to the SDP bound, which is more tight, they need more time to do so, and they also need more time to achieve the best Θ (see the columns "time to achieve Θ "). On the unconstrained quadratic problems and heaviest subgraph problems, our solver clearly outperforms SB and CSDP: in table 4, the gain in time is by a factor of 3 compared to SB and a factor of 2 compared to CSDP. Note that the gain is less important for graph bisection problems: the average time is still better, but CSDP performs better for some graphs (ex84f in particular). We highlight the quick convergence of our solver by plotting the run for a characteristic instance of each problem: in Figure 2 for a quadratic problem, in Figure 3 for a bisection problem, and similar plots for heaviest subgraph problems are in [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF] and [START_REF] Malick | Numerical study of semidefinite bounds for the k-cluster problem[END_REF].

The figures illustrate that our solver has a very strong initial increase effect. In particular, the improvement of the bound is very large in the first iterations. Moreover, we can observe in the tables that the running times of our solver are similar for instances of the same size. These two properties are highly desirable in the context of branch-and-bound. Each of the two SDP solvers has one of the two properties mentioned above. CSDP is also quite robust with respect to computing time (in contrast with SB for which "accidents" occurred, see for example ex204a in table 6). On the other hand, SB also shows a nice initial increase, although the computing times of SB are in general larger than those of CSDP. We notice in figure 3 that the initial heuristic of SB for finding a starting point is very efficient.

We point out that during the solving process, current primal and dual solutions of CSDP are not feasible and thus the dual value is not necessarily a valid lower bound. Nevertheless, we have observed during all our experiments that CSDP gives in fact lower bounds. We see also on Figure 3 that CSDP deteriorates the dual bound in order to make a big improvement in the feasibility of the current dual solution in the middle of the run.

Conclusions, perspectives

This paper introduces a family of semidefinite bounds for 0-1 quadratic programming. It presents a methodology to construct these bounds and to use them in practice. Numerical experiments illustrate the properties of these bounds on three combinatorial optimization problems.

The approach could be applied to any 0-1 problem of the form (1), but it is certainly not pertinent in any situation. Branch-and-bound algorithms using linear programming have reached a high level of sophistication and performance; they can solve efficiently many problems of the type (1). They have nevertheless shown limitations on some problems (such as the three problems considered here), and consequently, alternative approaches have been recently developed using quadratic programming (in particular the so-called QCR method, see [START_REF] Billionnet | Using a mixed integer quadratic programming solver for the unconstrained quadratic 0-1 problem[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]) or semidefinite programming (see [START_REF] Rendl | Solving max-cut to optimality by intersecting semidefinite and polyedral relaxations[END_REF]).

Our approach is relevant in these latter situations. We have a parameter α to set the targeted level of tightness below the usual SDP bound. Experiments with α small show that we can (almost) attain SDP bounds quicker than state-of-the-art SDP solvers. In addition, preliminary computational experiments with α large show that we can attain linear programming bounds and quadratic programming bounds quicker than CPLEX. Although this computational study has yet to be completed, we can say that this new family of bounds have some universality: whatever type of bound one looks for (tight, cheap, or balanced), one may be able to find it in this family.

Tests such as those of Section 5.2 are useful for setting the control parameter α to the desired tightness level. Once α is set, the bounds are computed by a quasi-Newton algorithm which combines several of the advantages of the SDP solvers: it is reliable (like CSDP) and has a strong initial increase effect (like SB). Those are properties of paramount importance for this method to be used as bounding procedure within branch-and-bound algorithms. A version of these bounds is used with success in [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF] in a branch-and-bound method to solve heaviest k-subgraph problems to optimality.

Finally, we mention that it is useful to combine these bounds with good heuristics for computing feasible solutions. Though cheaper than the usual SDP bounds, computing the bounds still demands sophisticated numerical linear algebra and optimization techniques. The bounds give better pruning in the search tree when combined with good feasible solutions, allowing us to fully benefit from the computing time investment in the bound quality.

Fig. 1

 1 Fig. 1 Comparison (axis: value of the bound vs cputime), for different values of α, of the bounds Θ uqp (µ, α) for a problem (4) of size n = 150.

Fig. 2 Fig. 3

 23 Fig. 2 An instance of the unconstrained quadratic problem (4) of size n = 100: comparison of the maximization of Θ uqp with the runs of SB and CSDP computing the SDP relaxation.

Table 2

 2 Numerical results for the heaviest subgraph problem (8) (each entry is the mean of 15 results). The table reports the computing time of Θ heav and gap(%), the averaged relative difference between the computed bounds and the best known solution of (8).

	Problem	Θ bis with α = 10 -2	Θ bis with α = 5.10 -3	Θ bis with α = 10 -3	Θ bis with 10 -4
	graph	k	time (s)	gap(%)	time (s)	gap(%)	time (s)	gap(%)	time (s)	gap(%)
	ex84a	42	0.09"	5.26%	0.11"	4.31%	0.20"	3.50%	0.29"	3.37%
		63	0.08"	7.12%	0.08"	6.02%	0.13"	4.93%	0.22"	4.74%
		49	0.11"	5.83%	0.14"	4.99%	0.15"	4.16%	0.26"	4.02%
		46	0.10"	5.58%	0.13"	4.63%	0.17"	3.81%	0.28"	3.67%
	ex84f	42	0.18"	1.33%	0.18"	1.26%	0.28"	1.20%	0.47"	1.20%
		63	0.15"	1.99%	0.18"	1.90%	0.29"	1.81%	0.45"	1.80%
		49	0.18"	1.54%	0.17"	1.47%	0.27"	1.39%	0.54"	1.38%
		46	0.18"	1.43%	0.21"	1.36%	0.35"	1.28%	0.49"	1.27%
	ex108a	54	0.18"	5.13%	0.18"	4.17%	0.24"	3.45%	0.40"	3.29%
		81	0.17"	6.23%	0.18"	5.03%	0.23"	4.04%	0.53"	3.29%
		63	0.16"	5.39%	0.18"	4.39%	0.27"	3.65%	0.52"	3.83%
		59	0.17"	5.34%	0.19"	4.37%	0.26"	3.64%	0.40"	3.48%
	ex132a	66	0.28"	4.67%	0.33"	3.66%	0.47"	2.92%	0.71"	2.76%
		99	0.24"	5.44%	0.27"	4.21%	0.40"	3.26%	0.81"	3.05%
		77	0.32"	4.88%	0.33"	3.89%	0.36"	3.12%	0.83"	2.96%
		72	0.28"	4.92%	0.36"	3.96%	0.42"	3.21%	0.57"	3.05%
	ex132f	66	0.51"	1.27%	0.58"	1.14%	1.07"	1.07%	1.34"	1.07%
		99	0.39"	1.65%	0.43"	1.54%	0.63"	1.45%	1.25"	1.43%
		77	0.41"	1.19%	0.52"	1.10%	0.68"	1.03%	1.22"	1.02%
		72	0.47"	1.23%	0.50"	1.08%	0.84"	1.01%	1.12"	1.00%
	ex156a	78	0.62"	5.06%	0.66"	4.10%	0.83"	3.36%	0.82"	3.22%
		117	0.34"	6.04%	0.47"	4.84%	0.61"	3.89%	1.23"	3.64%
		91	0.40"	5.15%	0.41"	4.19%	0.61"	3.43%	0.97"	3.28%
		85	0.46"	4.79%	0.52"	3.85%	0.58"	3.11%	1.00"	2.96%
	ex180a	90	0.82"	4.72%	1.04"	3.78%	1.27"	3.00%	1.71"	2.83%
		135	0.46"	5.36%	0.53"	4.19%	0.74"	3.25%	1.52"	3.02%
		105	0.54"	4.79%	0.69"	3.79%	0.81"	3.01%	1.93"	2.84%
		98	0.62"	4.71%	0.86"	3.73%	1.17"	2.97%	2.20"	2.81%
	ex204a	102	1.16"	4.06%	1.50"	3.11%	1.58"	2.36%	2.07"	2.21%
		153	0.96"	5.32%	0.72"	4.13%	1.19"	3.21%	2.55"	3.03%
		119	0.90"	4.50%	1.10"	3.51%	1.36"	2.74%	2.46"	2.60%
		111	0.95"	4.10%	1.07"	3.12%	1.68"	2.36%	2.61"	2.23%
	mean		0.40"	4.28%	0.46"	3.49%	0.62"	2.87%	1.06"	2.67%

Table 4

 4 Numerical tests for the unconstrained quadratic problem (each entry is the mean of 15 results).

	Problem	Θ heav with α = 10 -4		SDP with SB		SDP with CSDP
	n	k	time (s)	gap(%)	time	time (s) to achieve Θ heav	time (s)	time (s) to achieve Θ heav
	80	20	0.18"	0.16%	0.93"	0.44"	0.33"	0.26"
		40	0.11"	0.05%	0.97"	0.66"	0.36"	0.30"
		60	0.15"	0.01%	1.38"	1.19"	0.34"	0.29"
	100	25	0.24"	0.17%	2.07"	0.84"	0.58"	0.45"
		50	0.16"	0.06%	1.98"	0.93"	0.56"	0.46"
		75	0.17"	0.01%	4.66"	3.12"	0.59"	0.50"
	300	75	3.62"	0.11%	21.49"	8.18"	9.73"	7.79"
		150	1.11"	0.04%	11.15"	10.13"	11.86"	9.52"
		225	3.05"	0.01%	34.02"	11.45"	8.78"	6.82"
	mean		0.97"	0.07%	8.74"	4.10"	3.68"	2.93"

Table 5

 5 Numerical tests for the heaviest subgraph problem (each entry is the mean of 15 results).

Table 6

 6 Numerical results for graph bisection problems.

Acknowledgement

We are grateful to Sourour Elloumi who provided us a lot of material for the numerical tests, and to Sofia Zaourar who helped us to develop parts of the solver for bisection problems.


Gwendoline Blanchet

Lionel Moisan

An explicit Sharpness Index related to Global Phase Coherence

Keywords: sharpness, global phase coherence, image restoration, blind deconvolution, image quality

come

INTRODUCTION

Deriving no-reference image sharpness metrics has several applications in Image Processing, in particular in image restoration (blind and non-blind deconvolution, removal of compression artifacts) and in image quality assessment (typically used to rate compression schemes). It was shown in [START_REF] Blanchet | Measuring the global phase coherence of an image[END_REF] that a sharpness metric, called Global Phase Coherence (GPC), could be defined by reference to the regularity (Total Variation) of random phase images. This metric was then used in [START_REF] Vu | S3: A spectral and spatial sharpness measure[END_REF] as a part of a more sophisticated sharpness measure, and in [START_REF] Wang | A joint image quality assessment method based on global phase coherence and structural similarity[END_REF] to improve the sensitivity of SSIM to blur and aliasing.

In practice, the usefulness of the GPC is limited by the fact that its computation requires several hundreds of Monte-Carlo simulations -each one requiring one Discrete Fourier Transform (DFT) -to estimate the first two moments of the distribution of the Total Variation of random phase images. In the present work, we propose a variant of the GPC obtained by considering, instead of random phase images, the equivalent Gaussian random fields, for which we can derive explicit formulas (Section 2). In Section 3, we define the Sharpness Index and compare it to GPC, before we present an application to non-parametric blind deconvolution in Section 4.

This work has been supported by the French National Research Agency under grant ANR-09-BLAN-0029-01.

TOTAL VARIATION OF GAUSSIAN FIELDS

Let u : Ω → R be a grey-level M × N image defined on the rectangular domain Ω = {0, . . . , M -1} × {0, . . . , N -1}. We call standard white noise a random image W : Ω → R such that all W (x) are independent random variables distributed according to a the Gaussian distribution N (0, |Ω|

(that is, mean 0 and variance (M N ) -1 ). Note that this normalization ensures that E W 2 = 1, where E denotes the expectation and W 2 = x∈Ω W (x) 2 . Combining u and W with a periodic convolution, we obtain a new random image U = u ⋆ W defined on Ω by

U (x) = y∈Ω W (y) u(x -y), (1)
where u means the Ω-periodization of u (that is, u is Nperiodic along each coordinate). It is easy to prove that this definition ensures that U is a stationnary Gaussian random field with mean 0 Ω and covariance matrix

E U (x)U (x ′ ) = 1 |Ω| y∈Ω u(y) u(y + x ′ -x).
We now consider the notion of (periodic) Total Variation, that associates to an image v : Ω → R the value

TV(v) := ∂ x v 1 + ∂ y v 1 = x∈Ω |∂ x v(x)| + |∂ y v(x)|,
where ∂ x v(x, y) = v(x + 1, y)v(x, y) and ∂ y v(x, y) = v(x, y + 1)v(x, y). This is a regularity measure that has been extensively used in image processing since the seminal paper of Rudin, Osher and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF]. In order to establish the main theoretical result of this paper (Theorem 1 below), we need to introduce the gradient cross-correlation function defined by

∀z ∈ Ω, Γ(z) = y∈Ω ∇ u(y) • ∇ u(y + z) T ,
where ∇ u = (∂ x u, ∂ y u) T is the (periodic) gradient of u. Note that Γ(z) is a 2 × 2 matrix, that can be written under the form

Γ(z) = Γ xx (z) Γ xy (z) Γ yx (z) Γ yy (z) .

Now we can state

Theorem 1 The expectation and the variance of TV(U ) are

E TV(U ) = (α x + α y ) 2|Ω| π , (2)
Var TV(U ) = 2 π z∈Ω α 2 x • ω Γ xx (z) α 2 x + 2α x α y • ω Γ xy (z) α x α y + α 2 y • ω Γ yy (z) α 2 y , (3)
where α x = ∂ x u 2 , α y = ∂ y u 2 , Γ is the cross-correlation matrix of ∇u and ω is the function defined by

∀t ∈ [-1, 1], ω(t) = t arcsin t + 1 -t 2 -1. ( 4
)
Proof: We have by linearity ∂ x U = (∂ x u)⋆W , so that the discrete random field ∂ x U is a stationnary Gaussian field whose marginal distributions have zero mean and variance

E(∂ x U (x)) 2 = y∈Ω 1 |Ω| (∂ x u(x -y)) 2 = α 2 x |Ω| ,
and then a standard computation leads to

∀x ∈ Ω, E |∂ x U (x)| = 2α 2 x π|Ω| .
Using a similar reasonning on ∂ y U , we obtain that the expectation of

TV(U ) = x∈Ω |∂ x U (x)| + |∂ y U (x)
| is given by (2) as announced.

We now consider the variance of TV(U ). We have

E TV(U ) 2 = a,b∈{x,y} x,y∈Ω E |∂ a U (x)∂ b U (y)|,
and the stationnarity of ∇U entails (with z = yx)

E TV(U ) 2 = |Ω| a,b∈{x,y} z∈Ω E |∂ a U (0)∂ b U (z)|. (5)
Each term of this double sum can be written under the form E |XY |, where (X, Y ) is a zero-mean 2-dimensional Gaussian vector with covariance matrix

E X 2 E XY E XY E Y 2 .
For the term of ( 5) corresponding to a = x and b = y for example, we have

X = ∂ x U (0) and Y = ∂ y U (z), thus E XY = E x,y∈Ω ∂ x u(-x)∂ y u(z -y)W (x)W (y) = 1 |Ω| x∈Ω ∂ x u(x)∂ y u(z + x) = 1 |Ω| Γ xy (z)
and the covariance matrix of (X, Y ) is

1 |Ω| α 2 x Γ xy (z) Γ xy (z) α 2 y
, so that by Lemma 1 below we obtain that

E |XY | = 2α x α y π|Ω| • ω Γ xy (z) α x α y , with ω(t) = cos(arcsin t) + arcsin t • sin(arcsin t) = √ 1 -t 2 + t arcsin t = ω(t) + 1.
Combining all terms arising from (5), we finally obtain that

E TV(U ) 2 = 2 π a,b∈{x,y} z∈Ω α a α b • ω Γ ab (z) α a α b , (6)
where the contributions of (a, b) = (x, y) and (a, b) = (y, x) are equal since Γ yx (z) = Γ yx (-z). In the end, the announced result arises from Var TV(U ) = E TV(U ) 2 -(E TV(U )) 2 , which simply amounts to change ω into ω in [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF].

Lemma 1 let Z = (X, Y ) T be a Gaussian random vector with zero mean and covariance matrix

E ZZ T = p 2 pq sin θ pq sin θ q 2 .
Then, one has

E |XY | = 2|pq| π (cos θ + θ sin θ).
The proof results from standard integration calculus.

A NEW SHARPNESS INDEX

Definition and relation to GPC

Using Theorem 1, we are now in position of defining a new Sharpness Index, inspired from the GPC developped in [START_REF] Blanchet | Measuring the global phase coherence of an image[END_REF].

Definition 1

The Sharpness Index of an image u : Ω → R is

SI(u) = -log 10 Φ µ -TV(u) σ , (7)
where µ = E TV(U ) and σ 2 = Var TV(U ) are defined in Theorem 1, and Φ(x) = (2π) -1/2 +∞ x e -t 2 /2 dt is the tail of the Gauss distribution.

As in [START_REF] Blanchet | Measuring the global phase coherence of an image[END_REF], we replace the uncomputable probability P{T V (U ) ≤ T V (u)} by its Gaussian approximation based on the first two moments of T V (U ) (which is reasonable because T V (U ) is very near to be Gaussian in practice). Thus, the fundamental difference between GP C(u) and SI(u) relies in the random field U that is used: random phase images for GP C(u), a Gaussian random field for SI(u). As we shall see further, these two random fields lead to very similar properties, but the great advantage of the Gaussian random field is that it provides explicit formulas for SI(u), whereas the estimation of GP C(u) requires heavy Monte-Carlo simulations. In that sense, one can see the Sharpness Index we just defined as an explicit ersatz of the GPC.

Numerical computation

As for the GPC, we have to consider two potential issues in the estimation of the Sharpness Index: the fact that we assume periodic images (because we use a periodic TV and periodic convolutions), and the fact that actual images are quantized (on 256 grey levels, or more). To avoid undesirable effects, we thus apply the pre-processing mentioned in [START_REF] Blanchet | Measuring the global phase coherence of an image[END_REF] before we compute the Sharpness Index. This leads us to the more robust Sharpness Index

SI ′ (v) = SI(Q(per(v))),
where Q is a dequantization operator (a (1/2, 1/2) Fourierbased translation), and per(v) the periodic component of v (avoiding periodization artefacts) defined in [START_REF] Moisan | Periodic plus smooth image decomposition[END_REF]. Note that the computation of the Fourier Transform of u = Q(per(v)) only requires two DFTs. As concerns the computation of SI(u), it involves the computation of the gradient ∇u = (∂ x u, ∂ y u) T , the norms α x = ∂ x u 2 and α y = ∂ y u 2 , and the cross-correlation images Γ xx , Γ xy , Γ yy , which can be quickly computed in Fourier Domain thanks to the relation

Γ ab = ∂ a u * • ∂ b u.
Overall, the computation of SI ′ (v) requires the computation of 6 DFTs (2 to compute û from v, 1 to get u, and 3 to get the cross-correlation images). In practice, the computation of SI ′ on a 512 × 512 image takes about 200 ms on a recent laptop with a C implementation (as a comparison, computing the GPC with N = 1000 Monte-Carlo samples takes about 1 minute).

Comparison between Sharpness Index and GPC

The aim of the present paper is to show that GPC can be advantageously replaced by the Sharpness Index without loosing its properties. From a theoretical viewpoint, it is interesting to notice that random phase images (obtained by phase randomization in Fourier Domain) and the equivalent Gaussian images (obtained by convolution with a white noise) produce very similar images visually. In fact, they only differ in Fourier Domain by a mutilplicative Rayleigh noise [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF]. Hence, one can expect that the Sharpness Index and the GPC will essentially measure the same thing. Indeed, SI(u) and GP C(u) are the same combination of three numbers: T V (u), µ = E T V (U ), and σ 2 = Var T V (U ), U being a different random field in each case. Experimentally, we observed that the values of µ obtained for GP C(u) and SI(u) (that is, the average total variation of a random phase image and the equivalent Gaussian image) were very close, with typical relative differences below 1%. As concerns σ, however, there is a real difference between the two random fields, the value being systematically higher for Gaussian fields (the ratio seems to be around 7-8 in general). Hence, we expect SI(u) to be significantly lower than GP C(u), but as we shall see now it does not have important consequences in practice, since both quantities behave the same with respect to image changes. To compare the behavior of SI(u) and GP C(u) as functions of u, we compute as in [START_REF] Blanchet | Measuring the global phase coherence of an image[END_REF] the 2D blur-noise diagram, that is, the function

f (r, σ) = SI ′ (g r ⋆ u + σN ),
where g r is the 2-D Gaussian convolution kernel with standard deviation r, and N a white noise image with unit variance in each pixel. This blur-noise diagram is compared to the one obtained for GPC in [START_REF] Blanchet | Measuring the global phase coherence of an image[END_REF] on the same image (Lacornou), and the results are shown on Fig. 1. As we can see, the absolutes values of SI(u) and GP C(u) differ, but their relative behavior with respect to noise and blur is quite similar.

APPLICATION TO BLIND DECONVOLUTION

In [START_REF] Blanchet | Measuring the global phase coherence of an image[END_REF], an application of the GPC to parametric blind deconvolution was presented. It was shown that an appropriate value of the standard deviation r of a Gaussian kernel could be selected in a Wiener blind deconvolution scheme by maximizing the GPC of the deconvolved image. We applied the same protocol by replacing the GPC with the Sharpness Index, and we found identical results, that is, the same value of r. Rather than reproducing the identical curve we obtained (up to the scale change we discussed above), we propose here to illus-trate the interest of the Sharpness Index by considering a more challenging problem: non-parametric blind deconvolution.

Instead of formulating blind deconvolution as an inverse problem associated to an unknwown convolution kernel, we propose to use a more direct approach: given an image u, find the symmetric convolution kernel k such that SI(k ⋆ u) is maximal. Solving this optimization problem is not easy, because the function F : k → SI(k⋆u) is neither concave nor smooth (due to the presence of the TV operator). This is why we adopt a very simple (though probably not optimal in terms of computation time) strategy, called stochastic minimization: start with k = δ 0 (Dirac mass in (0,0)) repeat N times:

define k ′ from a random perturbation of k if SI(k ′ ⋆ u) > SI(k ⋆ u), then k ← k ′ return k ⋆ u
In order to limit the dimensionality of the problem, we considered separable kernels k with a small (21×21) support. Thus, considering the symmetry constraint on k, the exploration space has dimension 10 (we chose to maintain k(0) = 1, which is equivalent to the normalization x k(x) = 1 since SI(λu) = SI(u) for any λ = 0). The algorithm above could be trapped in a local maximum of the map F , but we did not observe this in practice. It might be that the map F associated to typical images only has one local maximum (which is global), but at this time we have no evidence of this. In the experiments we performed on several classical images, we used N = 10, 000, which seemed sufficient since we observed little improvement on the Sharpness Index after N = 1000. Also, running the algorithm several times always led to the same solution (up to a good level of precision).

The deconvolution algorithm we just presented was applied to several classical images, without applying any prior blur or noise. The results were surprisingly good, considering that no regularity constraint was imposed on the blur kernel or on the deconvolved image, as is generally done in classical blind deconvolution algorithms. Notice, however, that for some images the result was a bit too much deconvolved (in the sense that a careful examination of the result revealed the presence of ringing around edges). A representative example is shown on Fig. 2.

CONCLUSION

While the Global Phase Coherence defined in [START_REF] Blanchet | Measuring the global phase coherence of an image[END_REF] requires heavy Monte-Carlo simulations to be estimated, the Sharpness Index we proposed here is computed exactly with only 6 Discrete Fourier Transforms. Though, these two image quality indices exhibit very similar behaviors, in particular a sensitivity to image sharpness in a strong sense, that is, by opposition to several causes of the loss of details in an image (not only blur, but also noise, aliasing, ringing). Hence, considering the Sharpness Index instead of GPC opens interesting perspectives, both from a theoretical viewpoint (the explicit formulas will probably ease further analyses) and as concerns applications (in particular image restoration), for which computation time and estimation errors are no more a limitation.

Fig. 1 .

 1 Fig. 1. 2D blur-noise diagrams of Lacornou image, corresponding to the Sharpness Index (top) and the GPC (bottom).

Fig. 2 .

 2 Blind deconvolution of Lena image. Selecting the separable symmetric 21 × 21 kernel k that maximizes the Sharpness Index of k ⋆ u (u being the original Lena image, bottom left) results in a sharper image (bottom right) that reveal some details of the original image, while keeping noise and ringing at an acceptable level.