Integer Programming Models for the Median of a 0-1 String Set Under Levenshtein Distance

Median and centre strings are fundamental in such applications as classification [14, 26], trace reconstruction and matrix completion [10, 24], DNA/RNA sequence alignment [22], character and shapes recognition [8, 11]. Various solution methods were developed for metrics based on the Hamming distance: from heuristics and approximation algorithms [12, 16, 17, 34, 35], to exact methods based on Integer Linear Programming (ILP) [3, 4, 28]. However, the Hamming distance is not well-suited to models arising in computational biology and, to some extent, information theory, where the Levenshtein distance [6] is often more appropriate. The increasing complexity introduced by this metric explains why algorithms for this problem are, at the moment, less numerous and with comparatively poorer performance.

Several heuristics were devised for the Median String Problem, such as perturbation iterative refinement algorithms. Starting from an initial string [1], edits are done one at a time in sequence [27] or multiple at once [9], the latter achieving a faster convergence at the expense of solution quality. The procedure is interrupted when no further improvements can be found. Critical aspects lie in the definition of good initial strings and of operation ranking [29]. Alternatively, combination iterative refinement algorithms make use of subsequent averaging steps among pairwise input strings, up to finding an output string.

To our knowledge, only one ILP model has been investigated so far [20]. In that paper, a finite set of strings is given, each one weighted by the probability of its occurrence. The edit costs $a,d,s$ are arbitrary and independent of the position of the character. The goal is to minimize the total cost paid to transform the median string into each of the strings given, weighted by string occurrence probabilities, minus an incentive (equal to $s$) for any character saved in the output string: in other words, the median string is penalized proportionally to its length. Although ingenious, this formulation appears of little practical benefit: up to 100 000 seconds of CPU are in fact reported to compute the median of ten 4-characters strings of length 10.

In this paper, we revise the ILP model proposed in [20] and compare it with two new ILP formulations for the Median String Problem. Unlike [20], we do not consider different string occurrence probabilities, nor different operation costs $a,d,s$, because flat coefficients in the objective function generally make the problem harder due to symmetry.

More general cases and the Centre String Problem will be approached and discussed in a subsequent study. Also, we focus on the simpler case where $𝒜$ is binary: an extension to non-binary alphabets is formally straightforward, but its efficiency remains to be ascertained and will be the subject of future research.

Limiting to a binary alphabet, and in particular encoding $𝒜$ in $\{0,1\}$, has the advantage of a numerical interpretation of the output string that allows us to directly use its 0-1 components in the formulation, eliminating the need for auxiliary variables. In fact, even for binary $𝒜$, in [20] the median string components are generically integer and not 0-1, and for that reason their formulation requires activation constraints containing a “big M”.

Note also that the median string length cannot be predicted in a closed form from the problem data, but is an indeterminate of the problem and, as such, it is introduced in [20] as a single integer decision variable of the model. We adopt a different approach and, in our formulations, “digitalize” this variable by means of a limited amount of 0-1 variables. Another improvement we propose regards the upper bound to the median string length: the value proposed in [20], i.e., the sum of the string lengths in $\mathrm{\Lambda }$, can in fact be replaced by twice the maximum length of a string in $\mathrm{\Lambda }$ minus 1.

Our computational findings demonstrate the advantage obtained with the modifications described.

The Levenshtein distance between two strings $\lambda \in {𝒜}^n$, $\rho \in {𝒜}^m$ can be computed by a simple recursion [33], which amounts to finding the shortest path from $00$ to $nm$ in a directed acyclic graph $G=(V,A\cup D\cup S)$ where

and $A,D,S$ respectively stand for add (or across), delete (or down) and substitute (or slant) arcs. All the arcs are directed from left to right and/or from up to down: those in $A$ link $h0$ to $h1$, $h1$ to $h2$ and so on, up right to $hn$ for all $h$; those in $D$ link $0k$ to $1k$, $1k$ to $2k$ and so on, down to $mk$; finally, those in $S$ diagonally link $hk$ to $(h+1)(k+1)$ (see Figure 1). Summarizing:

We identify arcs with transitions that transform $\lambda$ into $\rho$ in the following way:

• $\mathrm{■}$

  An across arc $hk\to h(k+1)\in A$ corresponds to updating $\lambda$ by adding ${\rho }_{k+1}$ immediately after its $h^{th}$ element.

• $\mathrm{■}$

  A down arc $hk\to (h+1)k\in D$ corresponds to deleting the ${(h+1)}^{th}$ element of $\lambda$.

• $\mathrm{■}$

  A slant arc $hk\to (h+1)(k+1)\in S$ corresponds to replacing ${\lambda }_{h+1}$ with ${\rho }_{k+1}$.

When node $nm$ is eventually reached, $\lambda =\rho$.

Graph $G$ has $3nm+n+m$ arcs. Let $a_e$ ($d_e,s_e$) be the cost of adding to $\lambda$ (deleting from $\lambda$, substituting) an entry in the transition identified by arc $e$. Then the weight of a shortest $(00,nm)$-path equals the weighted Levenshtein distance between $\lambda$ and $\rho$. Defining $y_e\in \{0,1\}$ for any $e\in E=A\cup D\cup S$, this distance is obtained by solving

where $E_u^{+}$ ($E_u^{-}$) is the set of arcs entering (leaving) the generic node $u$, with $E_{00}^{+}=E_{nm}^{-}=\mathrm{\varnothing }$. Since the polyhedron is integral, $y_e$ does not need to be constrained in $\{0,1\}$.

The following definition is provided in [23] and [25] among many others.

Besides the method seen in § 2.1, we can represent the problem of finding the Levenshtein distance between $\lambda \in {𝒜}^n,\rho \in {𝒜}^m$ as that of finding a planar matching $M$ in a complete undirected bipartite graph $B=(L\cup R,L\times R)$ with node sets $L,R$ corresponding to the components of $\lambda ,\rho$, $|L|=n,|R|=m$. Precisely:

• $\mathrm{■}$

  Matching $h\in L$ to $k\in R$ corresponds to substituting ${\lambda }_h$ with ${\rho }_k$.

• $\mathrm{■}$

  An unmatched node of $R$ in position $k$ corresponds to adding to $\lambda$, in that position, the corresponding element of $\rho$.

• $\mathrm{■}$

  An unmatched node of $L$ corresponds instead to deleting the corresponding element from $\lambda$.

Since the matching is planar, the letters are inserted, replaced, or removed in the same order as they appear in the sequences. Then, to get the weight of the operations required to transform $\lambda$ into $\rho$, one simply has to add the weight of the unmatched nodes to that of the edges, which represent costly substitutions (that is, those with ${\lambda }_h\ne {\rho }_k$).

As an example, consider the planar matching depicted in Figure 2 and interpreted in the caption. This matching corresponds to the shortest path shown in Figure 1.

Defining $y_{hk},l_h,r_k\in \{0,1\}$ for all $h\in L,k\in R$, with

$y_{hk}=1$ if and only if $hk\in M$;

$l_h=1$ if and only if $h\in L$ is left unmatched;

$r_k=1$ if and only if $k\in R$ is left unmatched;

the distance between $\lambda$ and $\rho$ is computed by solving

where

• $\mathrm{■}$

  $s_{hk}$ is the cost of replacing ${\lambda }_h$ with ${\rho }_k$;

• $\mathrm{■}$

  $d_h$ ($a_k$) is the cost of deleting ${\lambda }_h$ from $\lambda$ (of adding ${\rho }_k$ to $\lambda$);

• $\mathrm{■}$

  contains all the edges that cross $uv$ in a canonical drawing of $B$.

The graph $B$ has $nm$ edges, as many as the essential variables of model (2) (the $l_h$ and $r_k$ variables are auxiliary and can be derived from the equations).

Let us begin with the distance computation presented in § 2.1. Let $n_i$ denote the length of ${\lambda }^i$ and suppose for the moment that we know the median string to be of length $m$. Associate then the $i^{th}$ string of $\mathrm{\Lambda }$ with a directed acyclic multigraph $G^i=(V^i,E^i)$ where

The only difference from the graph of § 2.1 is that each slant arc is duplicated: that is, we define $S^i=C^i\cup F^i$, where $C^i$ ($F^i$) is the set of costly (free) slant arcs of $G^i$. Consistent with our initial assumption of unweighted distance, we let $c_e=0$ for $e\in F^i$ and $c_e=1$ in all remaining cases. Given its shape, we call $G^i$ a grid (see Figure 1) with nodes arranged into rows (indexed by $h$) and columns (indexed by $k$).

The path used in § 2.1 to compute the distance between $𝐱$ and ${\lambda }^i$ is identified by a unit flow that traverses $G^i$ from source 00 to sink $n_{i}m$. Let $y_e\in I{R}_{+}$ denote the flow in arc $e$, and call $E_{hk}^{i+}$ (respectively, $E_{hk}^{i-}$) the set containing the arcs that enter (leave) node $hk$. Considering also duplicated slant arcs, $E_{hk}^{i+}$ and $E_{hk}^{i-}$ contain up to four arcs; in particular, $E_{00}^{i+}$ and $E_{n_{i}m}^{i-}$ are empty. Finally, if $e=(h-1)(k-1)\to hk\in F^i$, then we let $h_e=h,k_e=k$.

To reduce the length of a path, one can take advantage of free arcs, provided that they have non-zero capacity. But the capacity of a free arc ending in $hk$ depends on ${\lambda }_h^i$ and $x_k$, being 1 if $x_k={\lambda }_h^i$ (because in this case replacement has no cost and one can use the free arc) and 0 if $x_k\ne {\lambda }_h^i$ (for in this case, replacement has a cost and the free arc must be interdicted). Thus, if ${\lambda }_h^i=1$, then the capacity will be $x_k$, while if ${\lambda }_h^i=0$, it will be $1-x_k$. This consideration leads to the following mixed-integer linear model written for $i=1,\mathrm{\dots },t$:

Setting $E={\bigcup }_{i=1}^t{E}^i$ the objective reads

The path formulation thus obtained is essentially derived from the model proposed in [20]. Equations (3) enforce that each path starts in 00 and ends in $n_{i}m$, and define the paths together with flow conservation conditions (4). Inequalities (5), (6), that strengthen (c1) and (c2) of [20], allow a path to use free slant arcs provided that $x_k={\lambda }_h^i$. Note that the $y$-variables need not be constrained in $\{0,1\}$, since the constraint matrix is totally unimodular for any choice of the $x$-variables in $\{0,1\}$.

However, a crucial detail omitted so far in our path formulation is that the median (or centre) string length is not known, which is why in [20] a further integer variable is introduced to represent it. This variable, named $m$, is subject to the constraints ($b$) of [20] that involve a so-called “big M”: in such cases, the bounds obtained by linear relaxation are generally very loose, and the formulation performance is consequently expected to be very poor. To limit this inconvenience, we digitalize variable $m$ by introducing the following 0-1 decision variables for , where $p$ and $q$ are suitable bounds to $m$:

• $z_k=1$ if the median (or centre) length $m$ is at least $k$;

• $z_k=0$ otherwise.

Naturally,

as $m\ge k+1$ implies $m\ge k$. The generic grid $G^i$ is extended up to the $q^{th}$ column by setting

and the flow $y_e$ in each arc $e$ entering column $k\in (p,q]$ is capped by

By (10), then, no flow can reach a node $n_{i}k$ for $k$ larger than the median length.

Since $m$ is unknown, also the sink node is not determined and depends on the $z_k$. The second condition (3), that assumes the sink in node $n_{i}m$, is then replaced by

(where we set by definition $z_p=1,z_{q+1}=0$). In fact,

• $\mathrm{■}$

  if $m=k$, then $z_k=1,z_{k+1}=0$: the right-hand side of (11) equals 1 and node $n_{i}k$ is fixed as the total flow terminal;

• $\mathrm{■}$

  if , then $z_k=0,z_{k+1}=0$: the right-hand side of (11) equals 0 and node $n_{i}k$ is not reached by any flow because of inequalities (10). In this case, then, the sink belongs to a column that precedes the $k^{th}$ one;

• $\mathrm{■}$

  if $m>k$, then $z_k=1,z_{k+1}=1$: again, the right-hand side of (11) equals 0 and so the flow that possibly enters $n_{i}k$ also proceeds out of it as a whole, which implies that $n_{i}k$ belongs to a column that precedes the sink.

We can prove the following upper bound on $m$.

The matching formulation is constructed after the distance computation method presented in § 2.2. For each ${\lambda }^i\in \mathrm{\Lambda }$ we start from a complete undirected bipartite graph $B^i=(L^i\cup R,L^i\times R)$, where nodes in $L^i$ and in $R$ are ordered, so that the $k^{th}$ node of $L^i$ corresponds to ${\lambda }_k^i$ and the $k^{th}$ node of $R$ to $x_k$. Because we do not know the median string length $m$, we suppose $|R|=q\ge m$. The basic idea is to associate an edge joining the $h^{th}$ node of $L^i$ to the $k^{th}$ node of $R$ with the operation of substituting ${\lambda }_h^i$ with $x_k$. To distinguish between the case ${\lambda }_h^i\ne x_k$ and ${\lambda }_h^i=x_k$ we duplicate the edges of $B^i$:

• $\mathrm{■}$

  $C^i=L^i\times R=$ the set of edges of $B^i$ that represent a costly substitution;

• $\mathrm{■}$

  $F^i=L^i\times R=$ the set of edges of $B^i$ that represent a free substitution.

Leaving the $h^{th}$ node of $L^i$ (the $k^{th}$ node of $R$) unmatched corresponds to bearing the cost of adding ${\lambda }_h^i$ to (of deleting $x_k$ from) the median string. Unmatched nodes are represented by self-loops added to $B^i$, so we define

• $\mathrm{■}$

  $H^i=\{hh:h\in L^i\}=$ the set of edges with both extremes in $L^i$;

• $\mathrm{■}$

  $K=\{kk:k\in R\}=$ the set of edges with both extremes in $R$;

and form in this way the final edge set

For every node of $L^i\cup R$ we then define the set of all edges (including self-loops) that have that node as an extreme:

Let $p$ be a lower bound to $m$, and let us introduce our 0-1 decision variables:

$y_e=1$ if and only if $e\in E^i$ is selected ($1\le i\le t$);

$x_k=$ the $k^{th}$ element of the median string ($1\le k\le q$);

$z_k=1$ if and only if the median string length is at least $k$ ().

That is, the $y$-variables are as the $y,l,r$ in § 2.2, while the $x$- and the $z$-variables are as in § 3.1. For $e=hk$, we let $h_e=h,k_e=k$. The feasible region is then described by the following conditions written for $i=1,\mathrm{\dots },t$:

Equation (14) forbids the selection of any edge that touches $k\in R$ if $k$ is larger than the median string length $m$. Condition (15) rewrites the no-cross condition in § 2.2, with $X_f^i$ being the set of edges of $E^i$ that cross $f$ in a canonical drawing of $B^i$: if $f$ is selected, then the right-hand side is 0 and no edge crossing $f$ can in turn be selected; otherwise, the right-hand side is $n_i$ since this is the largest number of edges in a matching of $B^i$. This condition can be further strengthened by redefining $X_f^i$ as follows for $f=uv$:

Finally, conditions (16), (17), (18) are equivalent to (5), (6), (9) of § 3.1.

The objective is formally written as (8), but now $c_e=d$ or $c_e=a$ for self-loops, $c_e=0$ for $e\in F^i$, and $c_e=s$ for all the other edges.

Table 1 reports the computational results of the three formulations on the rnd dataset. The data reveal clear differences among the performance of the formulations as the maximum string length $n$ and the pool size $t$ increase. All in all, pat is consistently the most efficient and scalable model, mat shows moderate performance for small and medium instances but quickly deteriorates as $n$ grows, while hay becomes impractical already for moderately large values of $n$ and $t$.

For short strings ($n=5$), all formulations solve every instance to optimality. However, runtime differs substantially: pat remains below $0.20$ seconds even for $t=50$, mat stays below $0.70$ seconds, whereas hay already requires up to $7.13$ seconds. Increasing the pool size from $10$ to $50$ affects hay much more severely than the other two formulations.

When the maximum string length allowed increases to $n=10$, the performance gap between formulations increases significantly. pat still solves all instances in $\le 3.1$ seconds, while mat takes up to $57.56$ seconds for $t=50$. The hay formulation shows a steep increase in runtime, reaching $205.83$ seconds for $t=40$ and hitting the time limit for $t=50$, where it solves two instances of three, with a very large gap (124.48%) on the unsolved one.

For $n=15$, the differences become even more evident. pat continues to solve all instances to optimality, with runtimes up to $30.09$ seconds. In contrast, for $t\ge 10$ mat could only solve one instance (RND-15-20), reaching the time limit in all the remaining cases and showing average gaps between $14.38\%$ and $35.68\%$, with a maximum of $46.96\%$ on rnd-15-50. The hay formulation performs better than mat for small $t$, solving all instances up to $t=30$, but reaches the time limit for $t\ge 40$, with gaps exceeding $100\%$.

The longest strings tested ($n=20$) clearly illustrate the scalability limits of the formulations. pat remains the only model capable of solving all instances within the time limit, with a runtime of up to $318.32$ seconds. Both mat and hay hit the time limit for every instance with $t\ge 20$ (mat also for $t=10$), but the gaps of hay are substantially larger for $t=40$ and $t=50$: mat shows average gaps around $50\%$ and a maximum of $60\%$, while hay reaches $116.03\%$. Figures 3–4 provide an illustrative overview of the statistics discussed.

Table 2 reports the computational results of the three formulations on the ecc dataset. We recall that the ecc instances contain codewords with fixed lengths and strong algebraic structure. This distinction has a clear impact on the behaviour of the formulations. Across all experiments, pat remains the most efficient and reliable model. It solves every instance to optimality, with modest increases in runtime as the pool size grows. For both Hamming and Reed-Muller instances, the runtime remains below $30$ seconds for $t=50$, and the optimality gap is always zero. This confirms that pat scales very well even when the underlying strings come from highly structured code families.

The performance of mat is more subtle. For the smallest Hamming instances ($t=10$), it solves all cases to optimality, albeit with runtime significantly larger than pat and hay. As the pool size increases, mat reaches the time limit more frequently, especially for Hamming codes with $t\ge 20$. Nevertheless, even when it cannot close the gap completely, mat often keeps it reasonably small (e.g., $5.30\%$ for $t=20$ and $17.95\%$ for $t=30$, on average, with a peak of $21.81\%$) and provides, in general, meaningful bounds and feasible solutions. A similar pattern appears in the rm dataset: although mat hits the time limit for all pool sizes, the resulting gaps remain relatively stable (between $19.24\%$ and $25.45\%$, with a peak of $28.42\%$), indicating a performance that degrades gradually rather than sharply.

The hay formulation has a different trend. It performs very well for small and medium Hamming instances ($t\le 30$), solving them all to optimality in a runtime that remains quite manageable (up to $569.50$ seconds for $t=30$). However, its scalability deteriorates for larger pool sizes: for $t\ge 40$, hay reaches the time limit and returns large optimality gaps (on average, $45.05\%$ and $50.57\%$, peaking at $112.14\%$). In rm, hay again solves all instances up to $t=30$, but its performance drops for larger pools: it solves two out of three instances for $t=40$ and none for $t=50$, with non-negligible gaps (although, in this case, smaller on average than mat). A graphical comparison of the mean CPU time and gap is provided in Figure 5.