Hugo Gimbert
email: hugo.gimbert@labri.fr

Youssouf Oualhadj
email: youssouf.oualhadj@labri.fr

Soumya Paul
email: soumya@imsc.res.in

Computing Optimal Strategies for Markov Decision Processes with Parity and Positive-Average Conditions

We study Markov decision processes (one-player stochastic games) equipped with parity and positive-average conditions. In these games, the goal of the player is to maximize the probability that both the parity and the positive-average conditions are fulfilled. We show that the values of these games are computable in polynomial time. We also show that optimal strategies exist, require only finite memory and can be effectively computed.

Introduction

Infinite turn-based games are useful tools for modeling and analyzing open reactive systems. Such games have been extensively studied, both in the deterministic and the stochastic setting. In the deterministic setting, there are usually 2 players, player 0 or the system player and player 1 or the environment player who take turns in making moves in an arena, thus producing an infinite sequence of moves called a play. The objective is given in terms of a Borel subset Φ of the set of plays. The aim of player 0 is to play in such a manner that the resulting play always lies in Φ. The aim of player 1 is the opposite: to foil player 0's aim and to make the play settle down outside Φ.

A correct controller of the system ensures the system behaves correctly in a demonic environment. Synthesizing such a controller amounts to computing a winning strategy σ for player 0: no matter how player 1 plays, player 0 always wins if she plays according to σ.

Markov decision processes. In systems where hardware failures and other random events occur, the behavior of the environment is typically represented as a stochastic process [START_REF] Kwiatkowska | Stochastic model checking. In Formal Methods for the Design of Computer[END_REF][START_REF] Mateescu | Sliding-window abstraction for infinite markov chains[END_REF]. Markov decision processes (MDP's) have proven to be a powerful [START_REF] Vardi | Multi-objective model checking of markov decision processes[END_REF][START_REF] Baier | Problem and verification of Markov decision processes[END_REF] yet algorithmically tractable [START_REF] Courcoubetis | Markov decision processes and regular events[END_REF] tool. In Markov decision processes, the environment no longer is a strategic player but its moves are chosen randomly according to fixed transition probabilities that depend on the current state of the system. This setting naturally extends to 2 1 2 player games but in this paper we focus on MDP's i.e. 1 1 2 player games. Winning almost-surely. In Markov decision processes there are generally few chances that player 0 has a strategy to win for sure, player 0 rather tries to maximize the probability to fulfill her winning objective. The best strategies for player 0 are almost-surely winning strategies, which guarantee a win with probability 1. The existence of almost-surely winning strategies is the solution concept we are focusing on in this paper.

The parity and positive average objectives. The class of ω-regular objectives is a sub-class of Borel objectives that naturally arise from logical specifications of reactive systems. The parity objective is an objective that is complete for ω-regular objectives, in the sense that it can express any such objective [START_REF] Grädel | Automata, Logics and Infinite Games[END_REF]. In a parity game, with every vertex is associated an integer called a priority. Player 0 has to play in such a way that the maximum priority visited infinitely often is even. The parity objective enjoys the luxury of being positionally determined, that is, either of the players always has a winning strategy that is independent of the history of the play (see e.g. [START_REF] Zielonka | Infinite games on finitely coloured graphs with applications to automata on infinite trees[END_REF][START_REF] Chatterjee | Quantitative stochastic parity games[END_REF][START_REF] Zielonka | Perfect-information stochastic parity games[END_REF]). The interest for parity games also stems from the close relation between parity games and µ-calculus [START_REF] Emerson | Tree automata, mu-calculus and determinacy[END_REF].

When dealing with performance evaluation, one needs to go beyond ω-regular objectives. This is needed in particular to measure the average performances of the system along a play. In MDPs equipped with mean-payoff or positive average objectives, every move is associated with a real number called a reward and Player 0 seeks to maximize the average reward along the play.

The Positive Average objective states that Player 0 is declared the winner if the average value of the flow of rewards is strictly positive. This slightly differs from the mean-payoff objective where player 0 rather tries to maximize the expectation of the average value of the rewards.

The definition of positive average and mean-payoff objectives seem very close to each other, yet they are quite different when turning to applications. The positive average condition can be more suited to express certain constraints on the Quality of Service (QoS) of critical systems. For example, it may be more desirable that a car speed controller reacts to any emergency event in 100 ms almost-surely rather than in 1 ms with probability 99% and 9000 ms with probability 1%. This preference is easily expressible using the positive average objective (reacting in less than 150 ms is fulfilled with probability 100% in the first case and only 99% in the second case) whereas the mean-payoff objective is maximized in the second case.

To guarantee both correct behavior and maximal performances it is necessary to consider boolean combinations of parity and positive-average conditions.

Our results. In this work we consider Markov decision processes equipped with the conjunction of parity and positive average conditions. Our main result is that in such games, the values are computable in polynomial time. Moreover, there exist optimal strategies with finite memory and these strategies are effec-tively computable. Finally, our definition of positive average MDPs is robust, since these results hold both under the liminf and the limsup semantics.

Comparison with previous work. Apart from parity positive average games, there has been other attemps to define winning conditions that mix both parity conditions and conditions about the mean value of rewards.

In mean-payoff parity games [START_REF] Chatterjee | Mean-payoff parity games[END_REF] the objective of the player is to maximize her long run average reward and satisfy the parity condition: in case the parity condition holds, the player's income is the mean value of rewards, in the opposite case it is -∞. This objective is hard to rephrase in the stochastic setting because to be able to compute expectations, one would have to arbitrarily choose some punishing constants to replace the -∞ payoff with. Moreover, as is shown by the car speed controller example, the positive average and the mean-payoff conditions may induce different preference orders on the set of outcomes. For these two reasons, it is hard to compare our results with the results of [START_REF] Chatterjee | Mean-payoff parity games[END_REF].

In an Energy Parity game [START_REF] Chatterjee | Energy parity games[END_REF], the player not only should win the parity game but she should ensure at the same time that her cumulative reward never goes below 0. It seems that in the deterministic setting the energy and the positive average conditions are closely related: if the player wins the positive average game, then there exists an initial amount of energy such that the player wins the energy game as well. However this is no more true in the stochastic setting, where we do not know any precise link between energy parity games and parity and positive average games.

Priority mean-payoff games [START_REF] Gimbert | Deterministic priority mean-payoff games as limits of discounted games[END_REF] are another attempt to mix parity and meanpayoff games. Here, the payoff associated to a play is the average value of the stream of payoffs seen at those moments where the highest priority seen infinitely often was visited. These games have some nice properties: there always exists positional optimal strategies in priority mean-payoff games and moreover these games are closely linked to discounted games [START_REF] Gimbert | Deterministic priority mean-payoff games as limits of discounted games[END_REF][START_REF] Gimbert | Limits of multi-discounted Markov decision processes[END_REF]. However, priority meanpayoff games seem quite incomparable with parity and positive average games, both algorithmically and semantically.

Organisation of the paper. In section 2, we introduce basic notions and some results about MDPs together with some material specific to the reachability and tail winning conditions. In section 3, we introduce parity and positive average games. In section 4, we present our main result (Theorem 4 and Corollary 3) about the polynomial time computability of values of parity and positive average games. Moreover we prove that finite memory strategies are sufficient (Theorem 5). In section 5, we prove that using the lim inf semantic does not change the values of parity and positive average Markov decisions processes.

Preliminaries

In this section we introduce several basic notions that will be used throughout the paper: game arenas, strategies, games and values and, finally, optimal, positively and almost-surely winning strategies.

Games and strategies

An arena describes the set of moves that players are allowed to use when playing the game.

Definition 1 (Markov decision processes and plays). A Markov decision process or arena is a tuple M = (V, V 0 , V R , E, δ) where: (V, E) is a directed graph without deadlocks, the set of vertices V is partioned into two parts, Player 0 vertices V 0 and chance vertices V R . We equip M with the function δ : V R → D(V ) that assigns to each chance vertex a probability distribution over V , such that δ(u)(v) > 0 if and only if (u, v) ∈ E. We denote by vE the set {u ∈ V | (v, u) ∈ E} of successors of a vertex v.

A play in the arena M is a sequence of vertices

v 0 v 1 • • • v n such that for every 0 ≤ i < n, (v i , v i+1 ) ∈ E.
Definition 2 (Subarenas). A subgraph (V , E ) of (V, E) is a subarena if and only if each vertex in V 0 has a successor in V , and all the successors of a vertex in V R are in V . We usually denote A[V ] the subarena obtained by the set V . Definition 3 (Games and winning conditions). Let V be a set of vertices. A winning condition is a Borel subset Φ ⊆ V ω . An infinite play is winning for Player 0 if it belongs to Φ. A game is a tuple (M, Φ) made of an arena M with vertices V and a winning condition Φ.

A strategy of the player tells her how to play the game. Formally, Definition 4 (Strategy). A strategy σ for Player 0 is a function σ : V * V 0 → V such that, for every finite play

w = v 0 v 1 . . . v k ∈ V * V 0 , σ(w) ∈ v k E.
Strategies with finite memory are finitely representable strategies.

Definition 5 (Finite memory and positional strategies). A strategy σ is finite memory if there exists a finite set M , the memory of the strategy, an initial memory m I ∈ M , and functions g :

M × V → M and h : M × V → V such that if v 0 . . . v k ∈ V * V 0 is a history and m 0 . . . m k+1 is a sequence determined by m 0 = m I and m i+i = g(m i , v i ) then σ(v 0 . . . v k ) = h(m k+1 , v k ).
A strategy σ is said to be positional if M is a singleton. A positional strategy σ can be specified just by a function from the set of vertices of Player 0 to the set of all vertices, that is, σ :

V 0 → V such that for all v ∈ V 0 , σ(v) ∈ vE.
In a Markov decision process M, once we have fixed a strategy σ for Player 0 and an initial vertex v, this defines naturally a probability measure P σ v on V ω .

Values, optimal strategies and almost-surely winning strategies

For a winning condition Φ ⊆ V ω , the value of a vertex v ∈ V with respect to strategy σ is denoted: val σ (q) = P σ v (Φ) , intuitively this is the probability that Player 0 wins if the play is consistent with the strategy σ.

Definition 6 (Values and optimal strategies). The value of a vertex is defined as

val(v) = sup σ val σ (v)
If a strategy σ for player 0 is such that val(v) = val σ (v) then σ is said to be optimal.

There is an alternative notion of solution to a stochastic game, which is more qualitative.

Definition 7 (Almost-surely and positively winning strategies). We say that Player 0 wins almost-surely (resp. positively) from a vertex v if she has a strategy σ such that P σ v (Φ) = 1 (resp. P σ v (Φ) > 0). In this case, vertex v is said to be almost-surely (resp. positively) winning for Player 0. The set of almost-surely (resp. positively) winning vertices for Player 0 is denoted W =1 (resp. W >0 ) and called the almost-surely (resp. positively) winning region of Player 0.

A positively winning vertex is also positively losing. Definition 8. A Markov decision process is qualitatively determined if every vertex is either almost-surely winning or positively losing for Player 0.

Reachability games

The simplest class of games is the class of reachability games. In a reachability game, the goal of the player is to reach a set of target states T ⊆ V , in other words the winning condition is V * T V ω .

In reachability games, the set of positively and almost-surely winning vertices is especially easy to compute, using elementary fixpoint algorithms.

Definition 9 (Positive and almost-sure attractors). The positive attractor for player 0 to a subset T of V is the set of vertices from which the player can reach T with positive probability. It is denoted by Attr >0 (T, V ) and is formally defined as follows. Let f : 2 V → 2 V be the operator such that for any U ⊆ V ,

f (U ) = {v ∈ V | vE ∩ U = ∅} . Then Attr >0 (T, V ) is the least fixed point of f (T ).
The positive attractor of the chance player to a subset T of V is denoted by Attr c >0 (T, V ) and is defined similarly.

The almost-sure attractor for player 0 to T is the set of vertices from which she can reach T with probability 1. It is denoted by Attr =1 (T, V ) and is defined as

Attr =1 (T, V ) = (V \ Attr c >0 (V \ Y, V )) ∪ T , where Y = Attr >0 (T, V ).
Note that there is a positional strategy for player 0 to "attract" the play from any vertex in Attr >0 (T, V ) (resp. Attr =1 (T, V )) to T with positive probability (resp. probability 1). Properties of positive and almost-sure attractors are given in the following proposition.

Proposition 1. The positive (resp. almost-sure) attractor of Player 0 to T is exactly the set of positively (resp. almost-surely) winning vertices of Player 0 in the reachability game.

The complement of a positive attractor for Player 0 is a trap for her.

Definition 10 (Safe set). Similar to the attractor set, we define the safe set of the player as the largest subarena from which the player has a strategy to avoid reaching T . Formally, it is denoted Safe(T, V ) and obtained as follow. Let

f : 2 V → 2 V such that for any subset U of V f (U ) = {v ∈ V p | vE ⊆ U } ∪ {v ∈ V c | vE ∩ U = ∅} ∪ T .
Then Safe(T, V ) is the complement of the least fixed point of f (T ).

Tail Games

We say that a game is tail if the winning condition that equip this game is tail.

Definition 11 (Tail winning condition). Let Φ ∈ V ω a winning condition. Φ is tail if ∀u ∈ V * and ∀w ∈ V ω then uw ∈ Φ ⇐⇒ w ∈ Φ .
For Markov decision processes equipped with tail winning condition, the notions of values and qualitative solutions are tightly linked. This is ilustrated by Theorem 1 which will be used throughout the paper. Theorem 1. [START_REF] Gimbert | Solving Simple Stochastic Tail Games[END_REF] In any Markov decision process equipped with a tail winning condition, the three following assertions hold.

1. Tail games are qualitatively determined. winning vertices. 2. One can implement an optimal strategy with the same memory as needed by an almost-surely winning strategy.

The vertices with value 1 are exactly the almost-surely

The above theorem shows that one can focus only on the computation of almost-surely winning regions when considering tail winning conditions on Markov decision processes. Proposition 2 is a general result that shows how to solve any disjunction of tail winning conditions. Proposition 2. Let Φ 0 , • • • , Φ n be n tail winning conditions and M a Markov decision process. The almost-surely winning region for the game

n i=0 Φ i is given by the set Attr =1 ( n i=0 W =1 [Φ i ], V ).
Proof. Let us prove that Player 0 has an almost-surely winning strategy from the set Attr =1 ( n i=0 W =1 [Φ i ], V ). Player 0 plays as follows. First she applies her attractor strategy until she reaches one of the W =1 [Φ i ] then she applies her almost-surely winning strategy σ i associated with the game Φ i . This strategy is clearly almost-surely winning since Player 0 reaches one of the W =1 [Φ i ] with probability 1.

To see that Player 0 cannot win almost-surely outside Attr =1 (

n i=0 W =1 [Φ i ], V ), let W = Attr =1 ( n i=0 W =1 [Φ i ],
V ) and consider the largest sub-arena A in V \W . There, Player 0 has no strategy to win almost-surely for any Φ i , by qualitative determinacy (see Theorem 1, fact 3) we get that ∀v ∈ A, ∀σ, ∀0 ≤ i ≤ n,

P σ v (Φ i ) = 0. This implies ∀v ∈ A, ∀σ, P σ v ( n i=0 Φ i ) = 0
, what shows that every vertex in A has value 0. For any other vertex v not in W and not in A, the probability for a given strategy σ that the play reaches A is strictly greater than 0 otherwise it would imply that v ∈ W , which terminates the proof.

Corollary 1. Let Φ a tail winning condition. Assume that W =1 [Φ] can be computed in polynomial time, then there exists a polynomial time algorithm to compute the value of each vertex.

Proof. The first fact of Theorem 1 implies that if Player 0 cannot win almostsurely then she cannot win positively, hence to win, Player 0 has to reach her almost-surely winning region. Using the third fact of the same theorem we get that in tail games the value of vertex is exactly the value of the reachability game played on the same arena and where the target set is the almost-surely winning region. Assuming that the almost-surely winning region can be computed in polynomial time and using the fact the reachability games can be solved in polynomial time finish the proof.

Parity Games and Positive-Average Games

In this section, we define parity games, positive-average games, and recall known results about these games.

Parity games

In parity games, the winner of the play is determined by the set of priorities visited infinitely often during the play.

Definition 12 (Parity games). Let C be a finite subset of N called the set of priorities. A parity game is played in a Markov decision process whose set of vertices V is labelled with a priority function χ : V → C that assigns to each vertex a priority. The parity winning condition is:

Par = {v 0 v 1 v 2 • • • ∈ V ω | lim sup n χ(v n ) is even} .
For any priority d ∈ C, we denote the set of vertices having priority d by

V d , that is, V d = {v ∈ V | χ(v) = d}.
Parity games are a very useful tool in verification and automata theory [START_REF] Grädel | Automata, Logics and Infinite Games[END_REF].

Theorem 2 ([23, 6, 8]). In Markov decision processes equipped with parity condition, Player 0 has a positional optimal strategy. Moreover, parity games are qualitatively determined and the almost-surely (resp. positively) winning regions are exactly the set of vertices with value 1 (resp. with strictly positive value). These values and the positional optimal strategies are computable in polynomial time.

Positive average games

In positive average games, Player 0 wants to maximize the probability that the average value of rewards is strictly positive.

Definition 13 (Positive average games).

A positive average game is played in a Markov decision process whose set of vertices V is labelled with a reward mapping r : V → Q that assigns to each vertex a rational number called the reward. The positive average winning condition is:

Avg >0 = v 0 v 1 v 2 • • • ∈ V ω | lim sup n→∞ 1 n n-1 i=0 r(v i ) > 0 . (1)
There is another natural definition of positive average games, which is very similar except the lim sup is replaced by lim inf. We denote this condition Avg >0 .

We shall show later (cf. Theorem 6) that the choice of either definition does not impact our results for one-player stochastic games.

Theorem 3. In a Markov decision process equipped with positive-average condition, Player 0 has a stationary optimal strategy. Moreover, parity games are qualitatively determined and the almost-surely (resp. positively) winning regions are exactly the set of vertices with value 1 (resp. with strictly positive value). These values and the positional optimal strategies are computable in polynomial time.

To prove this theorem we use mean-payoff games and use the fact that these games can be solved in polynomial time using linear programming. Definition 14 (Mean-payoff Games). A mean-payoff game is played in a Markov decision process whose set of vertices V is labelled with a reward mapping r : V → Q that assigns to each vertex a rational number called the reward. The value of a vertex v in a mean-payoff game is

val(v) = sup σ E σ v lim sup n→∞ 1 n n-1 i=0 r(v i ) .
Proof (Theorem 3). Since positive-average games are tail, the first part of the theorem follows directly from Theorem 1.

Since the set of almost-surely winning vertices is exactly the set of vertices of value 1, computing the values amounts to compute the almost-surely winning region. In order to compute this set we first use linear programming to compute the set of vertices W with value > 0 for the mean-payoff game. This can be done in polynomial time [START_REF] Putterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]. Let τ an optimal strategy for the mean-payoff game, by [START_REF] Liggett | Stochastic games with perfect information and time average payoffs[END_REF][START_REF] Gillette | Stochastic games with zero stop probability[END_REF] we know that this strategy can be chosen memoryless, let also M[τ ] the Markov chain induced by τ . In M[τ ] consider C the set of all recurrent states. We show that Player 0 wins the positive-average game almost-surely from any vertex in W = Attr =1 (C, V ). Let π the attraction strategy and let σ the following strategy, Let us show that Player 0 cannot win almost-surely outside W . We know that from any vertex outside W , the probability to reach a recurrent state c such that E τ c lim sup n→∞

σ(v) = π(v) if v ∈ W \ C and σ(v) = τ (v) if v ∈ C.
1 n n-1
i=0 r(c i ) ≤ 0 is strictly greater than 0, hence the value of every vertex outside W is striclty less than 1. Corollary 1 concludes the proof.

Corollary 2. In any Markov decision process M we have:

∀v ∈ V, val Avg >0 (v) = val Avg >0 (v) .
The proof of this corollary is postponed to Section 5 where it is proved for a larger class of games.

4 Solving Parity and Positive-average Games with lim sup semantics

Computing the Values

In this section we consider Markov decision processes equipped with Par ∧ Avg >0 winning condition. We give a polynomial time algorithm that computes the almost-surely winning region and we show that Player 0 needs memory of size exponential in the size of the arena. We characterize the winning regions by induction on the priorities available in the arena. The two following lemmata proved in the appendix characterize the almost-surely winning regions when the highest priority is even (Lemma 1) and when the highest priority is odd (Lemma 2). Lemma 1. In any Markov decision process M where the winning condition is Par ∧ Avg >0 and the highest priority d is even, the set of almost-surely winning vertices is given by the largest set W ⊆ V such that The core idea of Lemma 1 is that in the set A[W \Attr >0 (V d ∩W, W )] Player 0 has an almost-surely winning strategy to achieve the Par ∧ Avg >0 objective, hence by being consistent with it she wins. In the set Attr >0 (V d ∩W, W )], Player 0 applies in turn the attraction strategy to ensure the parity objective and her almost-surely winning strategy for Avg >0 objective until her accumulated reward gets strictly greater than some integer. By doing so Player 0 is ensured that after one alternation of the attraction strategy and the Avg >0 strategy, the sum of her average reward is always strictly positive. Let us show that Algorithm 1 is correct. First it considers the largest subarena where the highest priority is 2, namely A[S]. In A[S], once a fixed point is reached, R satisfies the conditions of Lemma 1. Finally it uses Lemma 2 to compute the almost-surely winning region in the whole arena. Using the fact that the original game can be rewritten as the disjunction of all the Φ d Proposition 2 shows that the almost-surely winning region is given by Attr =1 ( d∈D W =1 [Φ d ]), where D is the set of even priorities.

1. A[W ] is a subarena of M.
We now argue on the running time complexity. Each Φ d can be solved in polynomial time, hence our procedure runs in O(|D| • L) where L is the time one needs to solve each Φ d .

From the above theorem and by Corollary 1, we get the following corollary.

Finite Memory strategies are sufficient

In this subsection, we examine the memory needed by Player 0 to win almostsurely. We prove the following theorem.

Theorem 5. For Player 0, memory of size exponential in the size of the arena is sufficient and necessary to implement her almost-surely winning strategies.

To prove this theorem, we first prove Proposition 3 where it is shown that almost-surely winning strategies with finite memory exist, in Proposition 4 we show that exponential size memory is sufficient and in the example of Figure 1 it is shown that exponential size memory is necessary.

The main idea behind the proof of Theorem 5 is that instead of applying her positive-average strategy until her reward goes above some well chosen integer, Player 0 will apply this strategy until she is ensured that her expected average reward becomes greater than some value.

Existence of Finite Memory Strategies

To establish the existence of finite memory strategies, we need the following two lemmata that show the relationship between the total-reward defined as follow lim inf n→∞ n i=0 r(v i ) and the average reward (see Definition 13).

Lemma 3 ([2]

). Let M be a finite irreducible Markov chain with reward. Let s be a recurrent state of M. Assume lim inf n→∞ n i=0 r(v i ) = ∞. Then there exists an η > 0 such that the average reward accumulated between two consecutive visits of s is at least η. Moreover η has polynomial bit complexity and depends on the number of states in M.

Lemma 4 ([2]

). In any Markov decision process M, Player 0 has a positional strategy σ such that

∀v ∈ W, P σ v lim inf n n i=0 r(v i ) = ∞ = 1 , ( 2
)
where W is the almost-surely winning region for the Avg >0 game.

Proof (Lemma 4). The winning condition Avg >0 is submixing and tail, hence there exists a positional optimal strategy [START_REF] Gimbert | Pure stationary optimal strategies in Markov decision processes[END_REF]. Therefore, there exists a positional almost-surely winning strategy. Thus by Corollary 2, σ is almost-surely winning for Avg >0 as well. Hence the following equation holds,

∀v ∈ W, P σ v lim inf n 1 n + 1 n i=0 r(v i ) > 0 = 1 .
The same strategy σ yields (2). Proposition 3. In any Markov decision process M equipped with Par ∧ Avg >0 objective, Player 0 has an almost-surely winning strategy with finite memory.

The main idea beyond Proposition 3 is that the strategy described in Section 4 can be decomposed into 2 phases.

Phase 1 when Player 0 is applying her attraction strategy π. Phase 2 when Player 0 is applying her Avg >0 strategy τ .

After these two phases, Player 0 is ensured that her average payoff is strictly positive.

To prove Proposition 3, we show that by considering only the expected average reward, one can bound the time that Player 0 has to play each of the strategies of Phase 1 and Phase 2 and yet fulfill almost-surely the Par ∧ Avg >0 objective.

Size of the Memory. The following proposition shows that an exponential memory is necessary and sufficient. Proposition 4. In a one-player Markov decision process M with the Par ∧ Avg >0 objective, memory of size exponential in the size of the arena is both necessary and sufficient to implement an almost-surely winning strategies.

Consider the Markov Decision Process depicted in Figure 1. We show that memory of size exponential in the size of the arena is necessary.

Player 0 wins almost-surely from any vertex v ∈ M. Let σ be an almostsurely winning strategy with finite-memory of size k. Denote T R the absorption time in vertex n in the Markov chain obtained by removing from M the edge (n, 0). E σ v [T R ] gives the expected time to reach n if the initial state is v. Thus the expected average reward for Player 0 on the path from vertex 0 to vertex n

is -E σ v [T R ]. We show by contradiction that k ≥ E σ v [T R ].
Since σ is almost-surely winning, the strategy cannot cycle forever in vertex n thus it leaves vertex n after at most k loops. Thus the expected accumulated reward on a path from vertex 0 to itself is less than k

-E σ v [T R ].
If this value was negative then according to the law of large numbers, the expected average reward would be almostsurely negative as well, a contradiction. Since E σ 0 [T R ] ≥ 2 n+1 1 -1 2 n (proof in appendix), then σ has a memory at least exponential in the size of the arena.

Combining results of Proposition 3 and Proposition 4 lead Theorem 5 and thus the following corollary. Corollary 4. In any Markov decision process M where the winning condition is Par ∧ Avg >0 , Player 0 has an optimal strategy with memory of size exponential in the size of the arena.

Solving Parity and Positive-average Games with lim inf semantics

In the previous section we studied parity and positive-average games with lim sup objective. An alternative definition of these games is to replace lim sup by lim inf.

We show that all results of the previous section hold for this alternative definition.

Definition 15.

Avg >0 = v 0 v 1 v 2 • • • ∈ V * | lim inf n→∞ 1 n n-1 i=0 r(v i ) > 0 .
To compute the values in Par ∧Avg >0 games, we use the previous result on optimality using finite memory and known results on Markov chains theory. Actually we show that the value of every vertex in the game Par ∧Avg >0 is equal to the one in the game Par ∧ Avg >0 . Proposition 5. In any Markov decision process M we have:

∀v ∈ V, val Par ∧ Avg >0 (v) = val Par ∧Avg >0 (v) .
To prove this Proposition 5 we need the following lemma: Lemma 5 (see e.g. [START_REF] Norris | Markov chains[END_REF]). Let M be a finite Markov chain and r : V → R a reward function. The following equality holds for almost all runs.

lim inf n→∞ n-1 i=0 r(v i ) n = lim sup n→∞ n-1 i=0 r(v i ) n .
Proof (Proposition 5). We show that the following inequalities holds:

∀v ∈ V, val Par ∧Avg >0 (v) ≤ val Par ∧ Avg >0 (v) . ( 3
) ∀v ∈ V, val Par ∧ Avg >0 (v) ≤ val Par ∧Avg >0 (v) . (4)
That ( 3) holds is trivial. It is a consequence of the fact that every winning strategy for Par ∧Avg >0 is also winning for Par ∧ Avg >0 .

To prove (4), notice that according to Corollary 4 Player 0 can play optimally using finite memory in the Par ∧ Avg >0 game, thus there exists a strategy σ which is optimal and with finite memory. Hence:

val Par ∧ Avg >0 (v) = P σ v (Par ∧ Avg >0 ) = P σ v (Par ∧Avg >0 ) ≤ sup σ P σ v (Par ∧Avg >0 ) = val Par ∧Avg >0 (v) ,
where the first equality is by definition of the value and the second is by Lemma 5. Therefore (4) holds and Proposition 5 is proved.

Proposition 5 leads the following theorem.

Theorem 6. In any Markov decision process M where the winning condition is Par ∧Avg >0 , the values are computable in polynomial time. Moreover the optimal strategies can be implemented using finite memory.

Conclusion

In this paper we have considered Markov decision processes equipped with parity and positive-average winning conditions. We have shown that finite memory optimal strategies exist and that a memory exponential in the size of the arena is sufficient to implement these strategies. We have also given a polynomial time algorithm to compute the values of such games. Moreover, our algorithm solves parity and positive-average games under both the lim sup and lim inf semantics.

A Proof of Theorem 2

Proof. Since parity games are tail conditions, the first part of the theorem follows directly from Theorem 1.

Let M be a Markov decision process, To prove the second part of the theorem, notice that the parity condition can be written as disjoint union of winning condition where in each one, Player 0 wins if she satisfies the parity condition played in a parity game with three priorities. In other words, we define Φ d as the game played in M where every vertex with even priority d is transformed into a vertex with priority 2, every vertex with priority greater than d is transformed to a vertex with priority 3 and every vertex with priority less than d is transformed into a vertex with priority 1. The game Φ d can be solved in polynomial time since the number of priorities is fixed. Since Player 0 wins the original parity game if and only if she wins one of the games Φ d for some even priority d, we get by Proposition 2 that the almost-surely winning region for the original parity game

W = Attr =1 ( d∈D W =1 [φ d ])
where D is the set of all even priorities. Corollary 1 shows that the values can be computed in polynomial time and thus finishes the proof.

B Proof of Lemma 1

Proof. To prove the lemma we show the following: (i) Any set X ⊆ V satisfying 1, 2 and 3 is almost-surely winning. (ii) The almost surely winning region satisfies 1, 2 and 3. We start by proving (i). Let X ⊆ V , we exhibit an almost-surely winning strategy σ for Player 0 from any vertex in X. Strategy σ is defined as follows. If the play ever reaches X \ Attr >0 (V d ∩ X, X), then Player 0 forgets everything that happened up to now and starts playing an almost-surely winning strategy σ =1 in A[X \ Attr >0 (V d ∩ X, X)], which exists by [START_REF] Brázdil | One-counter stochastic games[END_REF]. As long as the play stays in Attr >0 (V d ∩ X, X), strategy σ is defined as follows. Initially, she plays her attractor strategy to V d ∩ X, π for |X| steps. Then she plays her almost-surely winning strategy for the Avg >0 objective, τ until her average payoff goes above an appropriately chosen η > 0. Then she starts from scratch.

We now prove that the strategy σ described above is almost-surely winning. If ever a play consistent with σ reaches X \ Attr >0 (V d , X), then according to Proposition 1 it will stay trapped in X \Attr >0 (V d , X). Thus any play consistent with σ is almost-surely winning.

Assume now that the play stays is in Attr >0 (V d , X). First, we show that the parity objective is satisfied. Let A n be the event: "A vertex with priority d is not visited within n back and forth switches between τ and π". Let q be the minimum of the probabilities of all the chance edges in the arena A. We have for every

v ∈ Attr >0 (V d , X) P σ v (A n ) ≤ (1 -q |X| ) n ≤ (1 -q |V | ) n . Now, because (1 -q |V | ) < 1, n>0 P σ v (A n ) ≤ n>0 (1 -q |V | ) n < ∞. Thus using Borel-Cantelli Lemma we get P σ v ( ∞ n=1 ∞ k=n A k ) = 0.
That is, the probability that infinitely many of the events A n happen is 0. Hence a vertex with priority d is eventually visited, and the parity objective is almost-surely satisfied.

Next, we prove that the positive-average objective is satisfied. Notice that since τ is positional, it also achieves the positive-average condition with the lim inf semantics. Thus there exists an integer η such that the accumulated average reward eventually never goes under η > 0. This is feasible since the play is happening in the almost-surely winning region for the positive-average condition. After reaching a vertex of priority d or applying the attractor strategy for |X| steps, Player 0 applies her strategy τ until her average payoff goes above η. Thus the Avg >0 objective is achieved almost-surely. The above facts show that σ is almost-surely winning. This completes the proof of (i). Note also that this strategy may require an infinite memory size since the time needed by Player 0 to make her average payoff greater than η is not bounded.

We now show (ii). Denote by W the almost-surely winning region. We prove that W satisfies 1, 2 and 3. That 1 holds is obvious. That W satisfies 2 follows from the fact that Player 0 can ensure almost-surely the conjunction Par ∧ Avg >0 in A[W ]. To see that 3 holds, note that G[W \ Attr >0 (V d , W )] is a trap for Player 0. So if she plays her almost-surely winning strategy σ defined on W , she wins almost surely the Par ∧ Avg >0 condition. This shows (ii).

C Proof of Lemma 2

Proof. We show that from any vertex in the described set, Player 0 has an almost-surely winning strategy. Player 0 applies the following strategy. As long as the game has not reached R, player 0 plays her attractor strategy π. If the play is in R, she uses her almost-surely winning strategy, τ , in R. That is,

σ :V -→ V σ(v) = π(v) if v ∈ R τ (v) if v ∈ R.
This strategy is almost-surely winning since any play consistent with it eventually reaches the set R and stays there forever. We now prove that the almost-surely winning region is exactly the set described by the lemma. For this, let W be a set of vertices satisfying the claim of the lemma. We show that Player 0 cannot win almost-surely in V \ W . Let σ a strategy and v ∈ V \ W . Any play consistent with σ either a) visits V \ Safe(V d , V ) infinitely often or b) the play ultimately reaches Safe(V d , V ). If case a) holds, using the Borel-Cantelli Lemma we get that a vertex of priority d is visited infinitely many times. If b) holds, let A be the random variable with values in |V | which denotes the first vertex reached in the set Safe(V d , V ). If P σ A (Par ∧ Avg >0 ) = 1, it would mean that A ∈ W . Hence P σ A (Par ∧ Avg >0 ) < 1 which implies P σ v (¬(Par ∧ Avg >0 )) ≥ P σ A (¬(Par ∧ Avg >0 )) > 0. This shows that the Player 0 cannot win almost-surely from any vertex in V \ W .

D Proof of Lemma 3

Proof. Let M be a finite irreducible Markov chain with reward. Suppose that the total reward of M diverge to infinity. Let s a recurrent state and T s the random variable that gives the time to the next visit of s when starting from s. We are interested in the following quantity E s 1 Ts Ts-1 i=0 r(v i ) . According to [START_REF] Brázdil | One-counter stochastic games[END_REF] (Lemma 11) this quantity is strictly positive if and only if the accumulated reward of M diverges to infinity. This proves the first part of the lemma.

We use a discounted approximation to compute a lower bound η. Let 0 < λ < 1 and V λ a vector such that for every state

s of M, V λ (s) = E s [(1 - λ) i≤0 λ i r(v i )]. It is well known that a) lim λ→1 V λ (s) = E s [ 1 n i≤0 r(v i )]
. This relation holds even in the more general case of simple stochastic games [START_REF] Gillette | Stochastic games with zero stop probability[END_REF][START_REF] Liggett | Stochastic games with perfect information and time average payoffs[END_REF]. b)V λ = (1 -λ)R + λP V λ , where R is the reward vector and P is the transition matrix of M. Hence, by b) we have that

V λ = (1 -λP ) -1 (1 -λ)R, thus V λ (s) = (1 -λP ) -1 (1 -λ)R(s).
This quantity is a rational fraction of λ, therefore there exists two polynomials Q and S with degree at most n, where n the number of states of M such that V λ (s) = Q(λ) S(λ) . According to a) we can write lim λ→1 V λ (s) = Q (1) S [START_REF] Baier | Problem and verification of Markov decision processes[END_REF] , this quantity has a polynomial-bit complexity. Thus, there is a polynomial P such that V λ ≥ 2 -P (n) . Using the strong Markov property we have lim λ→1

V λ (s) = E s [ 1 ts ts-1 i=0 r(v i )].

E Proof of Proposition 3

Proof. Let M be a Markov decision process. We prove by induction on the number of priorities that Player 0 has an almost-surely winning strategies with finite memory. Suppose that M has one priority c. If c is even then Player 0 plays a positive average game, according to Theorem 3, there exist a positional optimal strategy for her. If c is odd then Player 0 has no winning strategy.

Suppose that Player 0 can win almost-surely using finite memory in any Markov decision process which contains less than d priority. Let M be a Markov decision process with d priorities.

If d is odd, according to Lemma 2, to win Player 0 applies her attractor strategy until she reaches the almost-surely winning region for the game Par ∧ Avg >0 played in the subarena A[Safe(V d , V )]. Note that in this subarena the number of priorities is strictly less than d and thus she has a finite memory strategy. Since the attraction strategy is memoryless, Player 0 has a finite memory strategy to win almost-surely if the highest priority is odd.

If d is even, according to Lemma 1, either Player 0 is playing in the almostsurely winning region for Par ∧ Avg >0 in the subarena A[V \ Attr 0 >0 (V d , V )] or she is playing outside that arena. In the former case, by induction, Player 0 has a finite memory strategy to win and the proof is done. In the latter case she applies her attractor strategy π for a specified time, then she switches to her positive-average strategy τ to ensure this objective as well. In the remaining of this proof, we are going to show that the time Player 0 should apply τ so that the positive-average objective is achieved can be bounded.

While being consistent with τ , Player 0 eventually reaches a recurrent state in the Markov chain M[τ ], denote this state u. Since τ is almost-surely winning for the Avg >0 objective, Lemma 4 shows us that the total reward of M will diverge to ∞. Lemma 3 on the other hand, shows that the expected average reward accumulated between two consecutive visits of u is strictly positive. Let η be the lowest expected average reward accumulated between two consecutive hits of a recurrent state in M[τ ] by Lemma 3 we know that η > 0. We recall also that the strategy that we are constructing does not differ much from the one of Lemma 1, where Player 0 plays in turn the strategy π and τ . Denote n 0 ≤ n 1 ≤ • • • the moments when Player 0 starts from scratch. To win almost-surely Player 0 has to ensure the following invariant throughout the play.

∀s ∈ V , ∀i, E σ s 1 n i+1 -n i ni+1-1 k=ni r(v k ) ≥ x > 0 . (5)
We show that if (5) holds then

∀s ∈ V , E σ s lim sup n→∞ 1 n n-1 k=0 r(v k ) > 0 .
Let x > 0, then for every vertex v we have

∀i, E σ s 1 n i+1 -n i ni+1-1 k=ni r(v k ) ≥ x ⇒ ∀i, E σ s 1 n i ni-1 k=0 r(v k ) ≥ x ⇒ lim sup n→∞ E σ s 1 n n-1 k=0 r(v k ) ≥ x (6) ⇒ E σ s lim sup n→∞ 1 n n-1 k=0 r(v k ) ≥ x . (7)
Where the trasformation from ( 7) to ( 7) is by Fatou lemma. Hence by using the strong law of large numbers we get that the positive-average objective is ensured almost-surely. After applying π for |V | steps, Player 0 starts applying τ and increments the counter i each time the recurrent state R is visited. Let n the number of visits such that Equation 5 holds. For each even priority d we need the following memory

M d = V ×{0, 1, 2}×{0, • • • , |V |-1}×{0, • • • , n}. Let Update : V ×M d → M d be the update function such that, (u, v, b, i, j) =                          (v, 0, i, j + 1) if (b = 0) ∧ (j < |V | -1) ∧ (χ(u) = d) . (v, 1, i, j) if (b = 0) ∧ [(j = |V | -1) ∨ (χ(u) = d)] . (v, 1, i, j) if (b = 1) ∧ (u ∈ R) . (u, 2, 0, j) if (b = 1) ∧ (u ∈ R) . (v, 2, i + 1, j) if (b = 2) ∧ (u = v) ∧ (i < n) . (v, 2, i, j) if (b = 2) ∧ (u = v) ∧ (i < n) . (v, 0, i, 0) if (i = n) .
The strategy σ : V × M d → V consists in applying π the attractor strategy whenever b = 0 and applying τ the Avg >0 strategy whenever b = 0.

F Proof of Proposition 4

Proof. We start with proving that exponential memory is sufficient. As in the proof of Proposition 3, let M a Markov decision process and by M[τ ], M[π] denote the Markov chains induced by τ and π respectively. We define the following random variables, Note that if all the rewards in M are strictly positive, Player 0 plays only for the parity objective. Hence no memory is required. Now assume that theres exist non positive rewards in M. We want to compute an integer n such that Equation [START_REF] Chatterjee | Mean-payoff parity games[END_REF] in the proof of Proposition 3 holds.

T j+1 -T j T n E σ s   1 T j+1 -T j Tj+1-1 i=Tj r(v i ) T j     = E σ s T n -T 0 T n E σ s 1 T 1 -T 0 T1-1 i=T0 r(v i ) ( 10
)
≥ ηE σ s 1 - T 0 T n (11) ≥ η 1 - E σ s [T 0 ] n (12)
≥ η 1 - E σ s [T 0 -|V |] + |V | n
Where the transformation from [START_REF] Gillette | Stochastic games with zero stop probability[END_REF] to [START_REF] Gimbert | Pure stationary optimal strategies in Markov decision processes[END_REF] holds by Lemma 3 and from [START_REF] Gimbert | Pure stationary optimal strategies in Markov decision processes[END_REF] to [START_REF] Gimbert | Solving Simple Stochastic Tail Games[END_REF] because T n ≥ n. Hence,

E τ s n-1 j=0 c j T n ≥ η 1 - E σ s [T 0 -|V |] + |V | n . (13)
From ( 8), ( 9) and [START_REF] Gimbert | Deterministic priority mean-payoff games as limits of discounted games[END_REF] we get

E σn s 1 T n Tn-1 i=0 r(v i ) ≥ |V | n m + E σ s [T R ] n m + η 1 - E σ s [T R ] + |V | n
Let us find a value for n such that

m n (|V | + E σ s [T R ]) + η n (n -E σ s [T R ] + |V |) > 0 We find n > E σ s [T R ] + |V | - m η (|V | + E σ s [T R ]) ≥ E σ s [T R ] + |V | -m2 Q(|V |) (|V | + E σ s [T R ])
We compute an upper bound for E σ s [T R ]. Let P , the sub-stochastic matrix obtained from the transition matrix of M[τ ] by replacing every recurrent entry by 0. From any state v in M[τ ], the time to absorption is given by (I -P ) -1 (v). Hence Using the same arguments as in Lemma 3, we get that this quantity is exponential in a polynomial in the size of the arena.

The proof that exponential size memory is necessary is given after Proposition 4 except for the lower bound on absorption time. This lower bound is established in Proposition 6.

G Details of Computations of Figure 1 Proposition 6. In the game depicted on Figure 1, the expected absorption time from vertex 0 to vertex n is at least

2 n+1 1 - 1 2 n . 22
Proof. We compute a lower bound for E σ v [T R ]. We use a first step analysis to compute E σ 0 [T R ]. We know that for every state 0 ≤ i ≤ n -1

E σ i [T R ] = 1 + 1 2 E σ 0 [T R ] + 1 2 E σ i+1 [T R ] .
and for i = n

E σ n [T R ] = 0 .
Thus we get

E σ 0 [T R ] = 2 n n-1 i=0 1 2 i = 2 n+1 1 - 1 2 n .

  σ is clearly positional. Let us show that σ almost-surely winning, from every vertex v ∈ W Player 0 reaches C with probability 1. Second from for each vertex c ∈ C we can write E τ c lim sup n→∞ 1 n n-1 i=0 r(c i ) > 0, since c is recurrent using the strong law of large numbers we get that P τ c lim sup n→∞ 1 n n-1 i=0 r(c i ) > 0 = 1, which shows that σ is almost-surely winning.

2 .

 2 Player 0 wins almost surely the Avg >0 played in the subarena A[W ], 3. Player 0 wins almost surely the Par ∧ Avg >0 game played in the subarena A[W \ Attr >0 (V d ∩ W, W )].

Lemma 2 .Theorem 4 .

 24 In any Markov decision process M where the winning condition is Par ∧ Avg >0 and the highest priority d is odd, the set of almost-surely winning vertices is given by Attr =1 (R, V ), where R is the almost-surely winning region for the Par ∧ Avg >0 game played in the subarena A[Safe(V d , V )]. In any Markov decision process M where the winning condition is Par ∧ Avg >0 , the almost-surely winning region is computable in polynomial time.Proof. To compute the values in polynomial time we use similar technics as in the proof of Theorem 2. For each even priority d, we create a new game Φ d played on the same arena where each priority d is transformed into 2, each priority greater than d is transformed into 3 and all the priorities less than d are transformed into 1. Then to solve each of these games we use the procedure described in Algorithm 1.

Corollary 3 .Algorithm 1

 31 In any Markov decision process M where the winning condition is Par ∧ Avg >0 , the values are computable in polynomial time. Computes the almost-surely winning region for Player 0 in one of the game Φ d 1 R ← ∅ 2 S ← Safe(V3, V ) 3 repeat 4 In the arena A[S] compute R the almost-surely winning region for Player 0 for Avg >0 game 5 In the arena A[R ] compute R the positively losing region for the parity game for Player 0. 6 S ← Saf e(R , S) 7 until R = ∅ 8 return Attr=1(R , V )

- 1 . 5 , - 1 . 5 , - 1 . 5 , - 1 . 5 , - 1 +1 - 1 Fig. 1 .

 15151515111 Fig.1. Player 0 wins in M if she visits state 0 infinitely many times and she satisfies the Avg >0 objective.

-

  T R : with values in N, is the time to absorption in M[τ ]. -T n : with values in N, is the time needed to visit the first recurrent state reached n + 1 times.

  Where the inequality holds because T n ≥ n and min v∈V {r(v)} is negative. Hence Where the first inequality holds because T n ≥ n and min v∈V {r(v)} is negative. Hence

  E σ s [T R ] ≤ |V | + max v∈|V | {(I -P ) -1 (v)} .


Inken Wohlers
email: inken.wohlers@cwi.nl

Rumen Andonov
email: randonov@inria.fr

Gunnar W Klau
email: gunnar.klau@cwi.nl

Algorithm engineering for optimal alignment of protein structure distance matrices

Keywords: Protein structure distance matrix alignment, Algorithm engineering, Integer linear programming, Branch-and-cut, Preprocessing, dali

Protein structural alignment is an important problem in computational biology. In this paper, we present first successes on provably optimal pairwise alignment of protein inter-residue distance matrices, using the popular dali scoring function. We introduce the structural alignment problem formally, which enables us to express a variety of scoring functions used in previous work as special cases in a unified framework. Further, we propose the first mathematical model for computing optimal structural alignments based on dense inter-residue distance matrices. We therefore reformulate the problem as a special graph problem and give a tight integer linear programming model. We then present algorithm engineering techniques to handle the huge integer linear programs of real-life distance matrix alignment problems. Applying these techniques, we can compute provably optimal dali alignments for the very first time.

Introduction

Proteins are chains of amino acid residues that fold into complex three-dimensional (3D) structures. These structures largely determine the biological function of a protein. Structural alignment is the problem of detecting structural similarities between two proteins-similarities that indicate, for example, a similar biological function or a common evolutionary origin. Due to numerous applications in biology and medicine, it is an extremely important problem in computational biology. Unlike the related problem of sequence alignment, for which polynomial algorithms exist, structural alignment is computationally difficult. Current structural alignment algorithms fall in one of two categories: (i) they aim at achieving a low RMSD (root mean square deviation) of the superpositioned substructures for, ideally, a large number of aligned residues, or, (ii) they are based on aligning the inter-residue distances of the two proteins. Given the distances between all residues of a protein, its 3D structure can, except for chirality, be precisely reconstructed [START_REF] Havel | The theory and practice of distance geometry[END_REF]. Inter-residue distance matrices are thus well-suited and widely used representations to align protein structures, and used in numerous programs.

From an algorithmic perspective, structural alignment algorithms can be divided into heuristics [START_REF] Di Lena | Fast overlapping of protein contact maps by alignment of eigenvectors[END_REF][START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF][START_REF] Kawabata | Protein structure comparison using the Markov transition model of evolution[END_REF] or exact methods [START_REF] Andonov | Maximum contact map overlap revisited[END_REF][START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF][START_REF] Malod-Dognin | Maximum cliques in protein structure comparison[END_REF][START_REF] Wohlers | Towards optimal alignment of protein structure distance matrices[END_REF][START_REF] Xie | A reduction-based exact algorithm for the contact map overlap problem[END_REF]. While the first are typically fast, they do not provide a theoretical quality guarantee of their results, let alone for optimality. Exact methods compute provably optimal solutions, but may suffer from performance problems. Among the exact approaches, methods that maximize the number of short aligned distances, or the contact map overlap (CMO), have been shown to be the fastest [START_REF] Andonov | Maximum contact map overlap revisited[END_REF][START_REF] Di Lena | Fast overlapping of protein contact maps by alignment of eigenvectors[END_REF]. However, the weak point of CMO approaches is that the scoring function is too basic [START_REF] Wohlers | Towards optimal alignment of protein structure distance matrices[END_REF]. The other extreme is the sophisticated dali (distance-matrix alignment) scoring function [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF]. Targeting higher precision, it takes into consideration all inter-residue distances and assigns individual scores. Yet, because of the computational complexity of structural alignment and the huge size of the problem instances, it is extremely difficult to design an exact algorithm that uses this or similarly complex scoring functions.

In this paper, we present first successes on computing provably optimal alignments for complete inter-residue distance matrices using the dali scoring function. We decide to focus on dali, because it is a widely used structural alignment program with a public web service and the accuracy of its alignments is higher or among those of the top methods in various benchmarks [START_REF] Kolodny | Comprehensive evaluation of protein structure alignment methods: scoring by geometric measures[END_REF][START_REF] Mayr | Comparative analysis of protein structure alignments[END_REF]. We introduce the structural alignment problem formally in Section 2, enabling us to express a variety of scoring functions used in previous work as special cases in a unified framework. In Section 3 we propose the first mathematical model for computing optimal distance matrix alignments. We therefore reformulate the problem as a special graph problem and give a tight integer linear programming (ILP) model. We present algorithm engineering techniques to solve this model to optimality for practical instance sizes (Section 4). Finally, in Section 5, we discuss our computational results that include first provably optimal dali alignments for non-trivial structural similarities.

Problem Formulation

We now formally define the problem of aligning two protein structure distance matrices. Our general definition will allow us to express many related problems that have been addressed in previous work as special cases.

A distance matrix D of a protein of length n is a symmetric n × n matrix. Rows and columns correspond to individual residues of the protein, in the order they appear in the protein sequence. The entries D i,j , 1 ≤ i ≤ n, 1 ≤ j ≤ n, denote the distances between residues i and j. Many variants are possible, e.g., distances between the C α or C β atoms, or the minimum distance between any atoms of the residues.

An alignment of two distance matrices A and B of two proteins of length n A and n B is a matching of a subset of residues from protein A with a subset of residues from protein B that respects the sequential order of residues in their respective chain. More precisely, it is a pair of sequences of indices (I, J) with I = (i 1 , i 2 , . . . , i nL ) and J = (j 1 , j 2 , . . . , j nL ) satisfying

|I| = |J| = n L , 1 ≤ i 1 < i 2 < . . . < i nL ≤ n A ,
and 1 ≤ j 1 < j 2 < . . . < j nL ≤ n B , with the interpretation that residue i p is aligned to residue j p for p = 1, . . . , n L . Such an alignment (I, J) induces an alignment of pairs of inter-residue distances and we define its score as

S(I, J) := nL o=1 nL p=1 s(A io,ip , B jo,jp ) + S gap (I, J), (1)
where s : R ≥0 × R ≥0 → R gives the scoring term of aligning individual interresidue distances from protein A with distances from protein B and S gap (I, J) is a function penalizing the gaps of the alignment (I, J). We denote by A(A, B) the set of all possible alignments of two distance matrices A and B.

The problem consists now in finding the best alignment between the two matrices with respect to S. It is NP-hard, following from [START_REF] Lathrop | The protein threading problem with sequence amino acid interaction preferences is NP-complete[END_REF], since variable length gaps are admitted into the alignment and interactions between amino acid residues from the sequence are admitted into the scoring function.

S(I, J) .

A large number of algorithms exists that address the distance matrix alignment problem using a variety of scoring functions, all of which can be expressed in our framework. In the following we present some of the most popular scoring schemes.

Contact map overlap (CMO).

Based on a distance threshold τ , two residues are defined to be in contact or not. Many algorithms have been presented [START_REF] Andonov | Maximum contact map overlap revisited[END_REF][START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF][START_REF] Xie | A reduction-based exact algorithm for the contact map overlap problem[END_REF] that aim at finding an alignment with the maximum number of shared contacts. This corresponds to S gap (I, J) = 0 and

s(A i,j , B k,l ) = 1 A i,j ≤ τ and B k,l ≤ τ 0 otherwise, for 1 ≤ i < j ≤ n A and 1 ≤ k < l ≤ n B and some τ ∈ R ≥0 .
dast scoring. Motivated by the fact that maximizing shared contacts often generates alignments with large RMSD, the dast algorithm [START_REF] Malod-Dognin | Maximum cliques in protein structure comparison[END_REF] focuses on local alignments with low RMSD. This corresponds to aligning cliques of distances with low mutual distance differences and can be expressed as S gap (I, J) = 0 and

s(A i,j , B k,l ) = 1 |A i,j -B k,l | ≤ τ -∞ otherwise, for 1 ≤ i < j ≤ n A and 1 ≤ k < l ≤ n B and some τ ∈ R ≥0 .
dali scoring. dali [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF] is a popular heuristic for protein structural alignment. Expressed in our model (1) its scoring function [START_REF] Hasegawa | Advances and pitfalls of protein structural alignment[END_REF][START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF] corresponds to S gap (I, J) = 0 and

s(A i,j , B k,l ) = 0.2 - |Ai,j -B k,l | 1 2 (Ai,j +B k,l ) e -( 1 2 (Ai,j +B k,l )/20) 2 i = j and k = l 0.2 otherwise.
In our experimental results we will focus on the dali scoring scheme.

matras scoring. Analogous to sequence substitution matrices, the structural alignment algorithm matras [START_REF] Kawabata | Protein structure comparison using the Markov transition model of evolution[END_REF], uses log-odds value matrices M as structural scores. A value M ⌊d1⌋,⌊d2⌋ indicates the log-likelihood that distance d 1 is aligned to distance d 2 . A positive log-likelihood means that the corresponding distances are aligned more likely than expected by chance. In our scoring model, matras scoring can be realized by using function S gap (I, J) for affine gap costs and s(A i,j , B k,l ) = M ⌊Ai,j ⌋,⌊B k,l ⌋ .

Protein threading problem (PTP). It is interesting to note that local PTP [START_REF] Collet | Local protein threading by mixed integer programming[END_REF] can also be modeled in this framework with similar log-odds-based scoring functions [START_REF] Marin | Frost: a filter-based fold recognition method[END_REF] as matras.

ssap scoring. ssap [START_REF] Taylor | Protein structure alignment[END_REF] is an early approach to structural alignment and uses a structural scoring function that can be expressed as

s(A i,j , B k,l ) = 50 (V A i,j -V B k,l ) 2 +2 i = j and k = l -(200 min{n A , n B }) 1 2
otherwise, where V A i,j and V B k,l are vectors between residue i and j in protein A and k and l in protein B, resp. Using the difference of vectors for scoring, ssap can account for the directionality of the compared inter-residue distances. Function S gap (I, J) assigns a penalty of -5 for each gap, independent of the gap's length.

Mathematical Model

We now present a mathematical model for structural alignment of protein distance matrices. To our knowledge, it is the first model for optimal protein structural alignment in a general setting as defined above. First, we reformulate the alignment problem in a graph-theoretical framework: An alignment graph G = (V, E) of two proteins of size n A and n B is a n A × n B grid graph in which each node i.k represents a possible alignment of a residue i in protein A to a residue k in protein B, see also Figure 1. Each row of G, from bottom to top, corresponds to a residue in protein A, and each column, from left to right, corresponds to a residue in protein B. Edges in G represent possible alignment of pairs of distances. An edge (i.k, j.l) denotes the possible alignment of the distance between residues i and j in protein A with the distance between k and l in protein B. Alignment graphs for structural alignment have first been used for CMO [START_REF] Andonov | Maximum contact map overlap revisited[END_REF], here we adapt them to a more general setting. We say a node j.l is strictly larger than i.k if and only if j > i and l > k and strictly smaller than i.k if and only if j < i and l < k. Edges are directed: any edge in the alignment graph has the property that its tail is strictly smaller than its head. It is thus oriented from south-west to north-east (cf. Figure 1). Because of the partial ordering between nodes and because there exist only edges between ordered nodes, the alignment graph is a directed acyclic graph. An alignment is a subset {i 1 .k 1 , i 2 .k 2 , . . .} of nodes that can be ordered such that each node is strictly larger than the previous one, i.e., i 1 < i 2 < . . . and k 1 < k 2 < . . .; we then say the nodes lie on a strictly increasing path. There is a one-to-one correspondence between alignments and strictly increasing paths. Every structural alignment can be represented by nodes on a strictly increasing path together with all induced edges (i.e. the associated augmented path). Its score is the sum of edge and node weights (cf. Figure 1, left). Determining the structural alignment of maximum score is equivalent to the maximum weight augmented path problem in the corresponding alignment graph.

α = 2s(A 1,2 , B 1,2 ), β = 2s(A 2,3 , B 2,4 ), γ = 2s(A 1,3 , B 1,4
Nodes in the alignment graph contradict if they are not on a strictly increasing path. A decreasing path contains a set of mutually contradicting nodes. That is, it is a set {i 1, right). We denote the set of all decreasing paths by C. Strictly smaller and strictly larger nodes exist only in strictly increasing paths and not in decreasing paths.

1 .k 1 , i 2 .k 2 , . . .} of nodes with i 1 ≥ i 2 ≥ . . . and k 1 ≤ k 2 ≤ . . . (cf. Figure
We introduce some further notation describing the neighborhoods of a node i.k. By V -(i.k) we denote the set of nodes that are strictly smaller than i.k (left neighborhood), and by V + (i.k) the set of nodes strictly larger than i.k (right neighborhood). The set C - i.k contains all decreasing paths in V -(i.k) and the set C + i.k all decreasing paths in V + (i.k). Our integer linear programming formulation uses two types of variables. Binary variables x ik represent nodes in the alignment graph and indicate whether residues i and k are aligned. Variables y ikjl denote whether the alignment graph edge (i.k, j.l) is present in the solution, i.e. whether distance A i,j is aligned with distance B k,l .

max nA-1 i=1 nA j=i+1 nB -1 k=1 nB l=k+1 2s(A i,j B k,l )y ikjl + nA i=1 nB k=1 s(A i,i B k,k )x ik (2) s.t. x ik ≥ j.l∈C y ikjl ∀C ∈ C + i.k , i ∈ [1, n A -1], k ∈ [1, n B -1] (3)
x ik ≥ j.l∈C

y jlik ∀C ∈ C - i.k , i ∈ [2, n A ], k ∈ [2, n B ] (4)
x ik ≤ 1 + j.l∈C s(Ai,j ,B k,l )≤0 (y ikjl -x jl ) ∀C ∈ C + i.k , i ∈ [1, n A -1], k ∈ [1, n B -1] (5) i.k∈C x ik ≤ 1 ∀C ∈ C (6)
x binary, y ≥ 0 [START_REF] Di Lena | Fast overlapping of protein contact maps by alignment of eigenvectors[END_REF] Note that the objective function (2) of this ILP only models the structural part of (1), and not S gap (I, J). It is possible to integrate also linear or affine gap costs. For the sake of simplicity, we do not describe this here but refer to [START_REF] Althaus | Multiple sequence alignment with arbitrary gap costs: computing an optimal solution using polyhedral combinatorics[END_REF] where this is done for multiple sequence alignment.

Constraints [START_REF] Collet | Local protein threading by mixed integer programming[END_REF] guarantee that the x variables form a proper alignment, i.e., a strictly increasing path in G. Constraints (3) and ( 4) link x and y variables. They prevent activating edges for which source or target node are not activated as well. Similarly, inequalities [START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF] force the activation of edges whose endpoints are activated. This is a novel class of constraints and necessary because edges (i.k, j.l) with s(A i,j , B k,l ) < 0 would otherwise never be part of an optimal solution. Note that all inequality classes have exponential size. Every feasible solution of model ( 2)-( 7) is a structural alignment and constitutes in the alignment graph a strictly increasing path with its induced edges.

Algorithm Engineering

In this section, we describe several algorithm engineering techniques to solve the integer linear program from the previous section to optimality for practically relevant instance sizes. We tried out several approaches, including solving a less tight but polynomially-sized model, developing a branch-and-cut approach for the full model, combinations of the two models, several preprocessing techniques, and a divide-and-conquer approach. Due to space limitations, we only describe the most successful approach here, which is a variable elimination preprocessing step followed by a branch-and-cut approach for the full model.

Preprocessing Using Variable Elimination

We reduce the number of alignment graph nodes by over-estimating the score of the best structural alignment including the respective node. If this overestimated score is less than the score of a known feasible solution, the corresponding node and its adjacent edges are fathomed, i.e., deleted from the alignment graph. The efficiency of this process depends on the quality of the used lower bound. We are currently using the heuristic solution provided by dali [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF].

In order to compute overestimations we use double dynamic programming, similar as in [START_REF] Andonov | Maximum contact map overlap revisited[END_REF]. In the first, local, dynamic program, we compute a profit p i.k for any node i.k ∈ V , which corresponds to the largest value that this node can add to the global score. To compute p i.k we focus only on incoming and outgoing edges from i.k, and consider only those with score s(A i,j , B k,l ) > 0. These scores are then assigned to the corresponding nodes j.l ∈ V -(i.k) and V + (i.k), respectively. The weight of the heaviest strictly increasing path

p + i.k in V + (i.k) (resp. p - i.k in V -(i.k)) can be computed in time proportional to the size of the rectangles V + (i.k) (resp. V -(i.k))
(i.k) V -(m.n) V + (m.n) .
Second from right: The same fine profit computation for filtering o.p. For filtering edge (m.n, o.p), we can take the minimum of the two cases. Right: Profit computation for i.k for fine filtering of edge (m.n, o.p).

The second level, global, dynamic program consists in overestimating any solution containing a given node m.n, i.e., in which x m.n = 1. For this purpose we associate profits p i.k to all nodes, and using them as weights we compute the heaviest strictly increasing path in V -(m.n) (respectively in V + (m.n)). Adding the weights of these two paths, as well as p m.n , we obtain an overestimated score for a node m.n.

The efficiency of this variable elimination depends on the preciseness of the profit computation. The coarsest but also fastest method to potentially eliminate a node m.n is to use the previously computed profits p i.k . The advantage is that the profits can be computed once for all nodes in the beginning. However, better upper bounds can be obtained by using profits p m.n i.k that are computed in the same way as p i.k but considering only nodes in V -(m.n) V + (m.n). This fine filtering gives tighter overestimations, but is more time-consuming.

Besides eliminating nodes, we can also eliminate edges. We implemented two ways of accomplishing this. One way, the so called coarse filtering of edge (m.n, o.p), uses the same profits as for coarse node filtering, cf. Figure 2 left. We can do better by re-using the profits p m.n i.k and p o.p i.k computed during the fine filtering for nodes (m.n) and (o.p). Instead of p i.k , we use in this case updated profits pi.k defined as pi.k = min(p m.n i.k , p o.p i.k ). If p m.n i.k and p o.p i.k are stored during the fine node filtering, such a finer edge filtering is then not more expensive than coarse filtering.

Branch-and-Cut

We use a cutting plane algorithm [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF] within a CPLEX branch-and-cut framework to solve ILP (2)- [START_REF] Di Lena | Fast overlapping of protein contact maps by alignment of eigenvectors[END_REF]. We first solve the LP relaxation of an initial ILP (with reduced number of constraints ( 3)-( 7)). In case of fractional solution, we solve a separation problem that yields a violated inequality cutting off this solution. In practice, often several of such cutting planes are generated. The new constraints are added to the relaxed ILP problem which is solved again, and so on. If no more cutting planes can be found, the formulation comprising the new constraints is solved by CPLEX within a branch-and-bound algorithm.

We now show how cutting planes for all types of constraints in model ( 2)-( 7) can be generated using Dijkstra's (or general shortest path) algorithm in a directed acyclic graph. For inequalities [START_REF] Collet | Local protein threading by mixed integer programming[END_REF] this is identical to [START_REF] Lenhof | A polyhedral approach to RNA sequence structure alignment[END_REF]. For the remaining inequalities the principle is the same.

Cutting planes for inequalities ( 6) can be generated efficiently by assigning value x ik obtained by the linear relaxation as weight to each alignment graph node. We can determine the most violated constraint by computing the maximum weight decreasing path C in the alignment graph nodes V . If this weight is greater than one, the nodes on the decreasing path violate constraint [START_REF] Collet | Local protein threading by mixed integer programming[END_REF]. Since node weights are between zero and one, we can always lift a constraint of type ( 6) by adding x-variables that refer to a corresponding decreasing path of maximum cardinality. Hence, we search for the maximum weight path amongst all decreasing paths of length (n A × n B ) -1. We reformulate this problem to a minimum weight path problem by assigning weights (1x ik ) to the nodes. Since all weights are positive, we can use Dijkstra's algorithm to solve this problem. If for the path C found in this way i.k∈C x ik > 1 holds, it generates a cutting plane.

Cutting planes of type ( 3) and ( 4) are generated analogous to cutting planes of type [START_REF] Collet | Local protein threading by mixed integer programming[END_REF]. For each node i.k we generate up to two cutting planes, one for incoming edges and one for outgoing edges. In the first case, we assign for fixed i.k the weights y jlik to each node j.l ∈ V -(i.k), and in the second case we assign y ikjl to each node j.l ∈ V + (i.k). Then we compute the maximum weight decreasing path in V -(i.k) and V + (i.k), respectively. If this path has weight greater than x ik , we identified a cutting plane. The reformulation of the maximum weight path problem to a minimum weight path problem that can be solved with Dijkstra's algorithm is analogous to cutting plane generation for constraints [START_REF] Collet | Local protein threading by mixed integer programming[END_REF]. In worst case we identify 2(n A -1)(n B -1) violated inequalities (3) and ( 4) in each iteration.

Similarly, we identify violated activation constraints [START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF]. For each node i.k, we assign the weight y ikjlx jl to each node j.l ∈ V + (i.k) with objective function coefficient s(A i,j , B k,l ) ≤ 0. Note that, because we consider only edges with negative objective function score, in this case we do not know the cardinality of the path beforehand. If the weight of the lightest decreasing path +1 is smaller than x ik , we identified a violated inequality [START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF]. Since we compute the minimum weight path in a directed acyclic graph with edge weights less than or equal to zero, we can not apply Dijkstra's algorithm. Instead we traverse all nodes in topological order, which is provided by sorting according to the above defined order on the nodes. A constraint of type (5) only cuts off the current solution if its x ik value is greater than zero. In practice, a good strategy is to generate cutting planes (5) only for nodes i.k with x ik = 1, because the alignment is mainly determined by edges (i.k, j.l) with positive score s(A i,j , B k,l ) and there are comparatively very few edges with negative score in a solution.

Results and Discussion

Data Sets and Setup of Computational Study

We evaluate variable elimination and the branch-and-cut approach on the Skolnick data set, a benchmark for clustering proteins. It consists of 40 protein chains of length between 97 and 255 residues that belong to five different SCOP (structural classification of proteins) families. This data set has been extensively used for the evaluation of CMO algorithms, e.g., in [START_REF] Andonov | Maximum contact map overlap revisited[END_REF][START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF][START_REF] Jain | Joining softassign and dynamic programming for the contact map overlap problem[END_REF][START_REF] Pelta | A simple and fast heuristic for protein structure comparison[END_REF]. dali skips hetero atom records and residues with incomplete backbone. Therefore we edited 7 PDB (protein data bank) files of the Skolnick data set; we changed heteroatom records to atom records and excluded residues that were for other reasons ignored during dali computation. Furthermore, in order to be in line with dali, we consider only the first decimal place of C α atom coordinates. The remaining computations were carried out with double precision leading to slight differences between the overall score reported by dali and the recomputed score.

To cluster the proteins in the Skolnick benchmark, we compute for each pair of proteins an alignment and a corresponding similarity score, resulting in 780 alignments. We compute them using dali and analyze the results. dali evaluates similarities using an empirical z-score which is based on the alignment's dali score. A z-score above 8 yields good structural superpositions [START_REF] Holm | Using Dali for structural comparison of proteins[END_REF]. In the Skolnick data set there are 164 alignments with z-score greater than 8. They all correspond to pairs of proteins from the same family and are considered "easy" instances in [START_REF] Andonov | Maximum contact map overlap revisited[END_REF]. All alignments between pairs of proteins from different families have a z-score of less than 4 (or dali detects no similarity). This gap between alignments with high z-score and alignments with low zscore is promoting the well-known fact that the Skolnick data set is a rather easy benchmark [START_REF] Di Lena | Fast overlapping of protein contact maps by alignment of eigenvectors[END_REF]. We focus in the subsequent analysis only on alignments with z-score greater than 8 because the poor performance of variable elimination for dissimilar proteins renders it currently infeasible to fit the problems into memory.

Furthermore, we use alignments from Sisy [START_REF] Andreeva | SISYPHUS-structural alignments for proteins with non-trivial relationships[END_REF][START_REF] Berbalk | Accuracy analysis of multiple structure alignments[END_REF], a more challenging data set that has been designed with the objective to provide difficult structural alignment instances. These alignments are individually inspected by experts and essentially manually created; therefore we consider them gold-standard reference alignments.

In order to compute optimal dali alignments using our branch-and-cut approach, we use CPLEX version 12.1 and a maximum running time of 30 hours. Instead of generating cutting planes for inequalities [START_REF] Collet | Local protein threading by mixed integer programming[END_REF], we use a polynomial number of constraints as initial ILP,

k l=1 x il + i-1 j=1 x jk ≤ 1, i ∈ [1, n A ], k ∈ [1, n B ] . (8)
These constraints from [START_REF] Andonov | Maximum contact map overlap revisited[END_REF] describe the set of increasing paths. We then add cutting planes (3), (4), and (5) within a CPLEX cutting plane callback function. Alignments have been computed on cluster nodes each equipped with two quad core 2.26 GHz Intel Xeon processors and 24 GB of main memory running 64 bit Linux.

Importance of Preprocessing

Real-life problem instances for aligning complete inter-residue distance matrices are huge, consisting of n A n B x-variables and nA 2 nB 2 y-variables. We therefore have to rely on a very effective preprocessing. In this work, we present variable elimination as preprocessing in order to minimize the number of variables that have to be considered explicitly in the dali ILP model. The effectiveness of this preprocessing depends on the similarity of the two proteins; if we apply preprocessing to two identical proteins, only the x-variables denoting identical residues and the y-variables of corresponding pairs of distances remain. As expected, the percentage of x-variables eliminated during preprocessing thus correlates with the dali z-score of the alignment (correlation coefficient 0.91, data not shown). The percentage of eliminated x-variables can be used successfully for correct classification of the Skolnick data set. This is an indication that highly similar proteins are easy to align optimally-an observation first made with exact CMO alignment algorithms.

Variable elimination has to be effective, firstly in order to reduce the model such that it fits into memory and secondly such that we do not have to add too many cutting planes. This is oftentimes feasible for the class of alignments that leads to good structural superposition (z-score greater than 8). For pairs of proteins with more subtle structural similarities, for example if only subsets of both proteins are structurally similar, this preprocessing fails. For these cases, the current variable elimination is a good first step, but it is not sufficient to reduce the model such that it can be handled by CPLEX.

Optimal Distance Matrix Alignments

Our branch-and-cut approach solves 75 of the 164 Skolnick alignments (46%) to provable optimality. The running times for the solved instances vary between 24s and 70296s and depend on the number of y-variables after variable elimination. The biggest solved instance has little less than 3 million y-variables. For 23 protein pairs (14%), CPLEX runs out of memory before the time limit is reached. Those instances have more than 15 million pairs of distances.

For 32 of the solved Skolnick instances (43%), the heuristic solution provided by dali was proven to be optimal. In the remaining instances, the optimal solution improved the heuristic solution slightly (less than 2% improvement in dali score). With improvements in this order of magnitude, dali might not fail to produce the optimal alignment, but determine a different alignment because it possibly uses different precision during computations. Furthermore, the compared Skolnick proteins are very similar in length and structure and the optimal alignment usually globally aligns them with a few small gaps. Therefore a good performance of the dali program is expected. We assume that in more difficult instances, in which the optimal alignment is not global, there is more room to improve on the heuristic dali solution.

We use small Sisy alignments in order to evaluate our branch-and-cut approach on such more difficult instances. Four instances are solved. We prove the optimality of the heuristic dali solution for two of them and for two improve with our exact solution the heuristic dali score slightly (by less than 1%).

For global alignments of structurally very similar proteins, e.g., Skolnick instances, the shift from a simple scoring function like CMO, to a more sophisticated scoring of all distances using for example the dali scoring function, shows little effect. Nonetheless, our previous evaluation of the complete Sisy data set illustrates that for detection of more complex structural similarities, basic scores as CMO are not always sufficient to obtain gold-standard reference alignments, and the dali scoring function performs on average better [START_REF] Wohlers | Towards optimal alignment of protein structure distance matrices[END_REF].

Conclusion and Future Directions

The branch-and-cut approach presented in this paper is the first exact algorithm that offers the feature of penalizing the alignment of non-compatible distances, e.g., a small and a large distance. Furthermore, it illustrates for the first time that it is feasible to compute provably optimal alignments of complete inter-residue distance matrices by means of preprocessing followed by a dedicated branch-and-cut approach that is implemented within a general purpose ILP solver. Especially, we were able to compute dali alignments to optimality, demonstrating that this popular and widely used heuristic structure comparison method generates optimal or close-to-optimal alignments, at least for relatively similar and relatively small problem instances.

The techniques we apply are currently successful only in the cases of similar proteins. Based on our experience, we believe that it is not realistic to provide exact solutions for all instances. On the other hand, all our results show that it is not necessary to consider all pairs of inter-residue distances in order to obtain good alignments. Methods from molecular distance geometry [START_REF] Wu | Rigid versus unique determination of protein structures with geometric buildup[END_REF], for example, can uniquely reconstruct a protein's 3D structure from a small subset of distances. Current alignments of sparse inter-residue distance matrices might thus perform promising because they capture to large extend the protein's structure. Nonetheless, we observe that it is essential to penalize non-compatible distances by assigning them a negative score. An illustration is given in Figure 3: the CMO scoring function greedily aligns as many residues as possible, which leads to an increase of alignment length at the expense of precise structural similarity of compact substructures as measured by RMSD.

Because of these observations, we plan to move towards structural alignment scoring schemes that reduce the size of the model by excluding pairs of distances. A first successful, yet still non-negative, scoring scheme was presented in [START_REF] Wohlers | Towards optimal alignment of protein structure distance matrices[END_REF]. Secondly, a purva, the implementation from [START_REF] Andonov | Maximum contact map overlap revisited[END_REF], is able to provide optimal CMO alignments even for very large and dissimilar instances. Based on these observations we will come up with hypotheses which pairs of distances can be excluded without losing alignment precision. Computing alignments to optimality, we can then objectively test these hypotheses.

Problem 1 (

 1 Distance matrix alignment) Given two distance matrices A and B, find an alignment (I, J) * of A and B with (I, J) * = arg max (I,J)∈A(A,B)

Fig. 1

 1 Fig.1Alignment graph for nodes that remain after preprocessing. Left: the set of nodes {1.1, 2.2, 3.4} forms a strictly increasing path that induces a structural alignment with score α + β + γ plus the score of the nodes. Here,α = 2s(A 1,2 , B 1,2 ), β = 2s(A 2,3 , B 2,4 ), γ = 2s(A 1,3 , B 1,4), and the node score s(A 11 , B 11 ) + s(A 2,2 , B 2,2 ) + s(A 3,3 , B 4,4 ). Right: Nodes on the decreasing path {4.1, 4.2, 3.2, 2.2, 2.3, 1.3, 1.4} mutually contradict.

Fig. 3

 3 Fig. 3 Alignment of 1aawA (gray) and 1gxiE (pink), an instance of the Sisy set. Optimal superposition according to the respective alignment. Residues colored in dark tone are aligned, residues colored in light tone are unaligned. Left: The Sisy reference alignment (29 aligned residues, RMSD of 1.14). Middle: The heuristic dali alignment correctly aligns all residues of the reference alignment, but extends the alignment length to 50 (RMSD of 2.55). Right: The optimal CMO alignment; it correctly aligns 96.55% of the aligned residues of the reference alignment. Alignment length is 56, RMSD 4.25. Additional gaps are inserted. Overaligning and insertion of additional gaps leads to a low RMSD value.

  using dynamic programming. Finally, we set p i.k = p - i.k + p + ik . Different computations for profit p i.k of node i.k. In black the feasible nodes. Left: Coarse profit computation. In black V -(i.k) and V + (i.k). Second from left: Fine profit computation for filtering node m.n when i.k < m.n. p - i.k stays unchanged, but p + i.k can be recomputed by considering only nodes in V +

			o		o
		m			m
	i	i	i		i
	k	k	n	k	p	k	n	p
	Fig. 2

Acknowledgements This work has been supported by DFG grant KL 1390/2-1 and was partly done when I. Wohlers was visiting IRISA supported by an INRIA grant. R. Andonov is supported by BioWIC ANR-08-SEGI-005 project and partially by DVU-01/197/16.12.2008 NSF-Bulgaria. Computational experiments were sponsored by the NCF for the use of supercomputer facilities, with financial support from NWO.