text:
# Improving Table Compression with Combinatorial Optimization An extended abstract appears in Proc. 13th ACM-SIAM Symp. on Discrete Algorithms, 2002.

## 1 Introduction

### 1.1 Table Compression

Table compression was introduced by Buchsbaum et al.  as a unique application of compression, based on several distinguishing characteristics. Tables are collections of fixed-length records and can grow to be terabytes in size. They are often generated by continuously operating sources and can contain much redundancy. An example is a data warehouse at AT&T that each month stores one billion records pertaining to voice phone activity. Each record is several hundred bytes long and contains information about endpoint exchanges, times and durations of calls, tariffs, etc.

The goals of table compression are to be fast, on-line, and effective: eventual compression ratios of 100:1 or better are desirable. While storage reduction is an obvious benefit, perhaps more important is the reduction in subsequent network bandwidth required for transmission. Tables of transaction activity, like phone calls and credit card usage, are typically stored once but then shipped repeatedly to different parts of an organization: for fraud detection, billing, operations support, etc.

Prior work  distinguishes tables from general databases. Tables are written once and read many times, while databases are subject to dynamic updates. Fields in table records are fixed length, and records tend to be homogeneous; database records often contain intermixed fixed- and variable-length fields. Finally, the goals of compression differ. Database compression stresses index preservation, the ability to retrieve an arbitrary record, under compression . Tables are typically not indexed at the level of individual records; rather, they are scanned in toto by downstream applications.

Consider each record in a table to be a row in a matrix. A naive method of table compression is to compress the string derived from scanning the table in row-major order. Buchsbaum et al.  observe experimentally that partitioning the table into contiguous intervals of columns and compressing each interval separately in this fashion can achieve significant compression improvement. The partition is generated by a one-time, off-line training procedure, and the resulting compression strategy is applied on-line to the table. In their application, tables are generated continuously, so off-line training time can be ignored. They also observe heuristically that certain rearrangements of the columns prior to partitioning further improve compression, by grouping dependent columns more closely.

We generalize the partitioning approach into a unified theory that explains both contiguous partitioning and column rearrangement. The theory applies to a set of variables with a given, abstract notion of combination and cost; table compression is a concrete case. To test the theory, we design new algorithms for contiguous partitioning, which speed training to work on-line on single files in addition to continuously generated tables; and for reordering in the off-line training paradigm, which improves the compression rates achieved from contiguous partitioning alone. Experimental results support these conclusions. Before summarizing the results, we motivate the theoretical insights by considering the relationship between entropy and compression.

### 1.2 Compressive Estimates of Entropy

Let $`𝒞`$ be a compression algorithm and $`𝒞(x)`$ its output on a string $`x`$. A large body of work in information theory establishes the existence of many optimal compression algorithms: i.e., algorithms such that $`|𝒞(x)|/|x|`$, the compression rate, approaches the entropy of the information source emitting $`x`$. These results are usually established via limit theorems, under some statistical assumptions about the information source. For instance, the LZ77 algorithm  is optimal for certain classes of sources, e.g., stationary and ergodic .

While entropy establishes a lower bound on compression rates, it is not straightforward to measure entropy itself. One empirical method inverts the relationship and estimates entropy by applying a provably good compressor to a sufficiently long, representative string. That is, the compression rate becomes a compressive estimate of entropy. These estimates themselves become benchmarks against which future compressors are measured. Another estimate is the empirical entropy of a string, which is based on the probability distribution of substrings of various lengths, without any statistical assumptions regarding the source emitting the string. Kosaraju and Manzini  exploit the synergy between empirical entropy and true entropy.

The contiguous partitioning approach to table compression  exemplifies the practical exploitation of compressive estimates. Each column of the table can be seen as being generated by a separate source. The contiguous partitioning scheme measures the benefit of a particular partition empirically, by compressing the table with respect to that partition and using the output size as a cost. Thus, the partitioning method uses a compressive estimate of the joint entropy among columns. Prior work  demonstrates the benefit of this approach.

### 1.3 Method and Results

We are thus motivated to study table compression in terms of compressive estimates of the joint entropy of random variables. In Section 2, we formalize and study two problems on partitioning sets of variables with abstract notions of combination and cost; joint entropy forms one example. This generalizes the approach of Buchsbaum et al. , who consider the contiguous case only and when applied to table compression. We develop idealized algorithms to solve these problems in the general setting. In Section 3, we apply these methods to table compression and derive two new algorithms for contiguous partitioning and one new algorithm for general partitioning with reordering of columns. The reordering algorithm demonstrates a link between general partitioning and the classical asymmetric traveling salesman problem. We assess algorithm performance experimentally in Section 4.

The new contiguous partitioning algorithms are meant to be fast; better in terms of compression than off-the-shelf compressors like gzip (LZ77); but not be as good as the optimal, contiguous partitioning algorithm. The increased training speed (compared to optimal, contiguous partitioning) makes the new algorithms usable in ad hoc settings, however, when training time must be factored into the overall time to compress. We therefore compare compression rates and speeds to those of gzip and optimal, contiguous partitioning. For files from various sources, we achieve 35–55% improvement in compression with less than a 1.7-factor slowdown, both compared to gzip. For files from genetic databases, which tend to be harder to compress, the compression improvement is 5–20%, with slowdown factors of 3–8.

The performance of the general partitioning with reordering algorithm is predicated on a theorized correlation between two measures of particular tours in graphs induced by the compression instances. We therefore measure this correlation, and the results suggest that the algorithm is nearly optimal (among partitioning algorithms). For several of our files, the algorithm yields compression improvements of at least 5% compared to optimal, contiguous partitioning without reordering, which itself improves over gzip by 20–50% for our files. In some cases, the additional improvement approaches 20%. While training time can be ignored in the off-line training paradigm, we show the additional time for reordering is not significant.

Finally, in Sections 57, we give some complexity results that link table compression to the classical shortest common superstring problem. We show that an orthogonal (column-major) variation of table compression is MAX-SNP hard when LZ77 is the underlying compressor. On the other hand, while we also show that the row-major problem is MAX-SNP hard when run length encoding (RLE) is the underlying compressor, we prove that the column-major variation for RLE is solvable in polynomial time. We conclude with open problems and directions in Section 8.

## 2 Partitions of Variables with Entropy-Like Functions

Let $`X=\{x_1,\mathrm{},x_n\}`$ be a set of discrete variables, each drawn from some domain $`𝒟`$, and consider some function $`H:𝒟^{}\mathrm{}`$. We use $`H(X,Y)`$ as a shorthand for $`H(Z)`$, where $`Z`$ is the set composed of all the elements in $`X`$ and $`Y`$: if $`X`$ and $`Y`$ are sets, then $`Z=XY`$; if $`X`$ and $`Y`$ are variables, then $`Z=\{X,Y\}`$; etc. For some partition $`𝒫`$ of $`X`$ into subsets, define $`(𝒫)=_{Y𝒫}H(Y)`$. We are interested in the relationship between $`H(X)`$ and $`(𝒫)`$. For example, let $`X`$ be a vector of random variables with joint probability distribution $`p(X)`$. Two vectors $`X`$ and $`Y`$ are statistically independent if and only if $`p(x,y)=p(x)p(y)`$, for all $`\{x,y\}`$; otherwise, $`X`$ and $`Y`$ are statistically dependent. Let $`H(X)=_{\{x_1,\mathrm{},x_n\}}p(x_1,\mathrm{},x_n)\mathrm{log}p(x_1,\mathrm{},x_n)`$ be the joint entropy of $`X`$. Then it is well known  that for any partition $`𝒫`$ of $`X`$, $`H(X)(𝒫)`$, with equality if and only if all the subsets in $`P`$ are mutually independent.

We can also view a table of $`n`$ columns as a system of $`n`$ variables. The relationship between certain compressors and entropy suggests that certain rearrangements that group functionally dependent columns will lead to better compression; Buchsbaum et al.  observe this in practice while restricting attention to partitions that preserve the original order of columns.

We are thus motivated to consider generally how to partition a system of variables optimally; i.e., to achieve a partition $`𝒫`$ of $`X`$ that minimizes $`(𝒫)`$, for some function $`H()`$, which we generally call the cost function. We introduce the following definitions. We call an element of $`𝒫`$, which is a subset of $`X`$, a class. We define two variables or sets of variables $`X`$ and $`X^{}`$ to be combinatorially dependent if $`H(X,X^{})<H(X)+H(X^{})`$; otherwise, $`X`$ and $`X^{}`$ are combinatorially independent. When $`H()`$ is the entropy function over random variables, combinatorial dependence becomes statistical dependence. Considering unordered sets implies that $`H(X,X^{})=H(X^{},X)`$. Note that in general it is possible that $`H(X,X^{})>H(X)+H(X^{})`$, although not when $`H()`$ is the entropy function over random variables. Finally, we define a class $`Y`$ to be contiguous if $`x_iY`$ and $`x_jY`$ for any $`i<j`$ implies that $`x_{i+1}Y`$ and a partition $`𝒫`$ to be contiguous if each $`Y𝒫`$ is contiguous. We now define two problems of finding optimal partitions of $`T`$.

###### Problem 2.1

Find a contiguous partition $`𝒫`$ of $`X`$ minimizing $`(𝒫)`$ among all such partitions.

###### Problem 2.2

Find a partition $`𝒫`$ of $`X`$ minimizing $`(𝒫)`$ among all partitions.

Clearly, a solution to Problem 2.2 is at least as good in terms of cost as one to Problem 2.1. Problem 2.1 has a simple, fast algorithmic solution, however. Problem 2.2, while seemingly intractable, has an algorithmic heuristic that seems to work well in practice.

Assume first that combinatorial dependence is an equivalence relation on $`X`$. This is not necessarily true in practice, but we study the idealized case to provide some intuition for handling real instances, when we cannot determine combinatorial dependence or even calculate the true cost function directly.

###### Lemma 2.3

If combinatorial dependence is an equivalence relation on $`X`$, then the partition $`𝒫`$ of $`X`$ into equivalence classes $`C_1,\mathrm{},C_k`$ solves Problem 2.2.

Proof.  Consider some partition $`𝒫^{}𝒫`$; we show that $`(𝒫)(𝒫^{})`$. Assume there exists a class $`C^{}𝒫^{}`$ such that $`C^{}C_i`$ for some $`1ik`$. Partition $`C^{}`$ into subclasses $`C_1^{},\mathrm{},C_{\mathrm{}}^{}`$ such that for each $`C_j^{}`$ there is some $`C_i`$ such that $`C_j^{}C_i`$. Let $`𝒫^{\prime \prime }=(𝒫^{}\{C^{}\})\{C_1^{},\mathrm{},C_{\mathrm{}}^{}\}`$. Since the $`C_i`$’s are equivalence classes, the $`C_j^{}`$’s are mutually independent, so $`H(C^{})_{j=1}^{\mathrm{}}H(C_j^{})`$, which implies $`(𝒫^{\prime \prime })(𝒫^{})`$. Set $`𝒫^{}𝒫^{\prime \prime }`$, and iterate until no such $`C^{}`$ exists in $`𝒫^{}`$.

If no such $`C^{}`$ exists in $`𝒫^{}`$, then either $`𝒫^{}=𝒫`$, and we are done, or else $`𝒫^{}`$ contains two classes $`C^{}`$ and $`D^{}`$ such that $`C^{}D^{}C_i`$ for some $`i`$. The elements in $`C^{}`$ and $`D^{}`$ are mutually dependent, so $`H(C^{},D^{})<H(C^{})+H(D^{})`$. Unite each such pair of classes until $`𝒫^{}=𝒫`$. $`\mathrm{}`$

Lemma 2.3 gives a simple algorithm for solving Problem 2.2 when combinatorial dependence is an equivalence relation that can be computed: partition $`X`$ according to the induced equivalence classes. When combinatorial dependence is not an equivalence relation, or when we can only calculate $`H()`$ heuristically, we seek other approaches.

### 2.1 Solutions Without Reordering

In the general case, irrespective of whether combinatorial dependence is an equivalence relation, we can solve Problem 2.1 by dynamic programming. Let $`E[i]`$ be the cost of an optimal, contiguous partition of variables $`x_1,\mathrm{},x_i`$. $`E[n]`$ is thus the cost of a solution to Problem 2.1. Define $`E[0]=0`$; then, for $`1in`$,

$$E[i]=\underset{0j<i}{\mathrm{min}}E[j]+H(x_{j+1},\mathrm{},x_i).$$
(1)
The actual partition with cost $`E[n]`$ can be maintained by standard dynamic programming backtracking.

If combinatorial dependence actually is an equivalence relation and all dependent variables appear contiguously in $`X`$, a simple greedy algorithm also solves the problem. Start with class $`C_1=\{x_1\}`$. In general, let $`i`$ be the index of the current class and $`j`$ be the index of the variable most recently added to $`C_i`$. While $`j<n`$, iterate as follows. If $`H(C_i\{x_{j+1}\})<H(C_i)+H(x_{j+1})`$, then set $`C_iC_i\{x_{j+1}\}`$; otherwise, start a new class, $`C_{i+1}=\{x_{j+1}\}`$. An alternative algorithm assigns, for $`1i<n`$, $`x_i`$ and $`x_{i+1}`$ to the same class if and only if $`H(x_i,x_{i+1})<H(x_i)+H(x_{i+1})`$. We call the resulting partition a greedy partition; formally, a greedy partition is one in which each class is a maximal, contiguous set of mutually dependent variables.

###### Lemma 2.4

If combinatorial dependence is an equivalence relation and all combinatorially dependent variables appear contiguously in $`X`$, then the greedy partition solves Problems 2.1 and 2.2.

Proof.  By assumption, the classes in a greedy partition correspond to the equivalence classes of $`X`$. Lemma 2.3 thus shows that the greedy partition solves Problem 2.2. Contiguity therefore implies it also solves Problem 2.1. $`\mathrm{}`$

### 2.2 Solutions with Reordering

Problem 2.2 asks for the best way to partition the variables in $`T`$, ignoring contiguity constraints. While a general solution seems intractable, we give a combinatorial approach that admits a practical heuristic.

Define a weighted, complete, undirected graph, $`G(X)`$, with a vertex for each $`x_iX`$; the weight of edge $`\{x_i,x_j\}`$ is $`w(x_i,x_j)=\mathrm{min}(H(x_i,x_j),H(x_i)+H(x_j))`$. Let $`P=(v_0,\mathrm{},v_{\mathrm{}})`$ be any path in $`G(X)`$. The weight of $`P`$ is $`w(P)=_{i=0}^\mathrm{}1w(v_i,v_{i+1})`$. We apply the cost function $`H()`$ to define the cost of $`P`$. Consider removing all edges $`\{u,v\}`$ from $`P`$ such that $`u`$ and $`v`$ are combinatorially independent. This leaves a set of disjoint paths, $`𝒮(P)=\{P_1,\mathrm{},P_k\}`$ for some $`k`$. We define the cost of $`P`$ to be $`(P)=_{i=1}^kH(P_i)`$, where $`P_i`$ is taken to be the unordered set of vertices in the corresponding subpath. If $`P`$ is a tour of $`G(X)`$, then $`𝒮(P)`$ corresponds to a partition of $`X`$.

We establish a relationship between the cost and weight of a tour $`P`$. Assume there are two distinct paths $`P_i=(u_0,\mathrm{},u_k)`$ and $`P_j=(v_0,\mathrm{},v_{\mathrm{}})`$ in $`𝒮(P)`$ such that $`u_k`$ and $`v_0`$ are combinatorially dependent and $`v_0`$ follows $`u_k`$ in $`P`$. In $`P`$ exist the edges $`\{u_k,x\}`$, $`\{y,v_0\}`$, and $`\{v_{\mathrm{}},z\}`$. We can transform $`P`$ into a new tour $`P^{}`$ that unites $`P_i`$ and $`P_j`$ by substituting for these three edges the new edges: $`\{u_k,v_0\}`$, $`\{v_{\mathrm{}},x\}`$, and $`\{y,z\}`$. We call this a path coalescing transformation. The following lemma shows that it is a restricted form of the standard traveling salesman 3-opt transformation, in that it always reduces the cost of a tour. It is restricted by the stipulation that $`u_k`$ and $`v_0`$ be combinatorially dependent.

###### Lemma 2.5

If $`P^{}`$ is formed from $`P`$ by a path coalescing transformation, then $`w(P^{})<w(P)`$.

Proof.  Consider

$$w(u_k,x)+w(y,v_0)+w(v_{\mathrm{}},z)$$
(2)
and

$$w(u_k,v_0)+w(v_{\mathrm{}},x)+w(y,z).$$
(3)
We have $`w(P^{})w(P)=(\text{3})(\text{2})`$. The definition of $`𝒮(P)`$ implies that $`(\text{2})=H(u_k)+H(x)+H(y)+H(v_0)+H(v_{\mathrm{}})+H(z)`$. That $`u_k`$ and $`v_0`$ are combinatorially dependent implies $`w(u_k,v_0)<H(u_k)+H(v_0)`$. Since $`w(X,Y)H(X)+H(Y)`$ for any $`X`$ and $`Y`$, we conclude that $`(\text{3})<(\text{2})`$. $`\mathrm{}`$

Repeated path coalescing groups combinatorially dependent variables. If a tour $`P`$ admits no path coalescing transformation, and if combinatorial dependence is an equivalence relation on $`X`$, then we can conclude that $`P`$ is optimal by Lemma 2.3. That is, $`𝒮(P)`$ corresponds to an optimal partition of $`X`$, which solves Problem 2.2. Furthermore, Lemma 2.5 implies that a minimum weight tour $`P`$ admits no path coalescing transformation.

When $`H()`$ is sub-additive, i.e., $`H(X,Y)H(X)+H(Y)`$, as is the entropy function, a sequence of path coalescing transformations yields a sequence of paths of non-increasing costs. That is, in Lemma 2.5, $`w(P^{})<w(P)`$ and $`(P^{})(P)`$. We explore this connection between the two functions below, when we do not assume that combinatorial dependence is an equivalence relation or even that $`H()`$ is sub-additive.

## 3 Partitions of Tables and Compression

We apply the results of Section 2 to table compression. Let $`T`$ be a table of $`n=|T|`$ columns and some fixed, arbitrary number of rows. Let $`T[i]`$ denote the $`i`$’th column of $`T`$. Given two tables $`T_1`$ and $`T_2`$, let $`T_1T_2`$ be the table formed by their juxtaposition. That is, $`T=T_1T_2`$ is defined so that $`T[i]=T_1[i]`$ for $`1i|T_1|`$ and $`T[i]=T_2[i|T_1|]`$ for $`|T_1|<i|T_1|+|T_2|`$. Any column is a one-column table, so $`T[i]T[j]`$ is the table formed by projecting the $`i`$’th and $`j`$’th columns of $`T`$; and so on. We use the shorthand $`T[i,j]`$ to represent the projection $`T[i]\mathrm{}T[j]`$ for some $`ji`$.

Fix a compressor $`𝒞`$: e.g., gzip, based on LZ77 ; compress, based on LZ78 ; or bzip, based on Burrows-Wheeler . Let $`H_𝒞(T)`$ be the size of the result of compressing table $`T`$ as a string in row-major order using $`𝒞`$. Let $`H_𝒞(T_1,T_2)=H_𝒞(T_1T_2)`$. $`H_𝒞()`$ is a cost function as discussed in Section 2, and the definitions of combinatorial dependence and independence apply to tables. In particular, two tables $`T_1`$ and $`T_2`$, which might be projections of columns from a common table $`T`$, are combinatorially dependent if $`H_𝒞(T_1,T_2)<H_𝒞(T_1)+H_𝒞(T_2)`$—if compressing them together is better than compressing them separately—and combinatorially independent otherwise.

Problems 2.1 and 2.2 now apply to compressing $`T`$. Problem 2.1 is to find a contiguous partition of $`T`$ into intervals of columns minimizing the overall cost of compressing each interval separately. Problem 2.2 is to find a partition of $`T`$, allowing columns to be reordered, minimizing the overall cost of compressing each interval separately. Buchsbaum et al.  address Problem 2.1 experimentally and leave Problem 2.2 open save for some heuristic observations.

A few major issues arise in this application. Combinatorial dependence is not necessarily an equivalence relation. It is not necessarily even symmetric, so we can no longer ignore the order of columns in a class. Also, $`H_𝒞()`$ need not be sub-additive. If $`𝒞`$ behaves according to entropy, however, then intuition suggests that our partitioning strategies will improve compression. Stated conversely, if $`H_𝒞(T)`$ is far from $`H(T)`$, the entropy of $`T`$, there should be some partition $`P`$ of $`T`$ so that $`H_𝒞(P)`$ approaches $`H(T)`$, which is a lower bound on $`H_𝒞(T)`$. We will present algorithms for solving these problems and experiments assessing their performance.

### 3.1 Algorithms for Table Compression without Rearrangement of Columns

The dynamic programming solution in Equation (1) finds an optimal, contiguous partition solving Problem 2.1. Buchsbaum et al.  demonstrate experimentally that it effectively improves compression results, and we will use their method as a benchmark.

The dynamic program, however, requires $`\mathrm{\Theta }(n^2)`$ steps, each applying $`𝒞`$ to an average of $`\mathrm{\Theta }(n)`$ columns, for a total of $`\mathrm{\Theta }(n^3)`$ column compressions. In the off-line training paradigm, this optimization time can be ignored. Faster algorithms, however, might allow some partitioning to be applied when compressing single, tabular files in addition to continuously generated tables.

The greedy algorithms from Section 2.1 apply directly in our framework. We denote by GREEDY the algorithm that grows class $`C_i`$ incrementally by comparing $`H_𝒞(C_iT[j+1])`$ and $`H_𝒞(C_i)+H_𝒞(T[j+1])`$. We denote by GREEDYT the algorithm that assigns $`T[i]`$ and $`T[i+1]`$ to the same class when $`H_𝒞(T[i,i+1])<H_𝒞(T_i)+H_𝒞(T[i+1])`$.

GREEDY performs $`2(n1)`$ compressions, each of $`\mathrm{\Theta }(n)`$ columns, for a total of $`\mathrm{\Theta }(n^2)`$ column compressions. GREEDYT performs $`2(n1)`$ compressions, each of one or two columns, for a total of $`\mathrm{\Theta }(n)`$ column compressions, asymptotically at least as fast as applying $`𝒞`$ to $`T`$ itself.

Even though combinatorial dependence is not an equivalence relation, we hypothesize that GREEDY and GREEDYT will produce partitions close in cost to the optimal contiguous partition produced by the dynamic program. We present experimental results testing this hypothesis in Section 4.

### 3.2 Algorithms for Table Compression with Rearrangement of Columns

We now consider Problem 2.2. Assuming that combinatorial dependence is not an equivalence relation, to the best of our knowledge, the only known algorithm to solve it exactly consists of generating all $`n!`$ column orderings and applying the dynamic program in Equation (1) to each. The relationship between compression and entropy, however, suggests that the approach in Section 2.2 can still be fruitfully applied.

Recall that in the idealized case, an optimal solution corresponds to a tour of $`G(T)`$ that admits no path coalescing transformation. Furthermore, such transformations always reduce the weight of such tours. The lack of symmetry in $`H_𝒞()`$ further suggests that order within classes is important: it no longer suffices to coalesce paths globally.

We therefore hypothesize a strong, positive correlation between tour weight and compression cost. This would imply that a traveling salesman (TSP) tour of $`G(T)`$ would yield an optimal or near-optimal partition of $`T`$. To test this hypothesis, we generate a set of tours of various weights, by iteratively applying standard optimizations (e.g., 3-opt, 4-opt). Each tour induces an ordering of the columns, which we optimally partition using the dynamic program. We present results of this experiment in Section 4.

## 4 Experiments

### 4.1 Data

We report experimental results on several data sets. The first three of the following are used by Buchsbaum et al. .

is a collection of 90-byte records from a customer care database of voice call activity.

is a collection of 32-byte records from a system of network status monitors.

is a portion of the United States 1990 Census of Population and Housing Summary Tape File 3A . We used field group 301, level 090, for all states. Each record is 932 byes.

is a file from Telcordia’s database describing local telephone switches. We appended spaces as necessary to pad each record to a uniform 30 bytes.

We also use several files from genetic databases, which are growing at a fast pace and pose unique challenges to compression systems . These files can be viewed as two-dimensional, alphanumeric tables representing multiple alignments of proteins (amino acid sequences) and genomic coding regions (DNA sequences).

The files EGF, LRR, PF00032, backPQQ, callagen, and cbs come from the Pfam database of multiple alignments of protein domains or conserved protein functions . Its main function is to store information that can be used to determine whether a new protein belongs to an existing domain or family. It contains more than 1800 protein families and has many mirror sites. The size of each table can range from a few columns and rows to hundreds of columns and thousands of rows. We have chosen multiple alignments of different sizes and representing protein domains with differing degrees of conservation: i.e., how close two members of a family are in terms of matching characters in the alignment.

The file cytoB is from the AMmtDB database of multi-aligned sequences of Vertebrate mitochondrial genes for coding proteins . It contains data from 888 different species and over 1100 multi-alignments of protein-coding genes. The tables corresponding to the alignments tend to have rows in the order of hundreds and columns in the order of thousands, much wider than the other files we consider. We have experimented with one multiple alignment: CytoB represents the coding region of the mitochondrial gene (from 500 different species) of cytocrome B.

Table 1 details the sizes of the files and how well gzip and the optimal partition via dynamic programming (using gzip as the underlying compressor) compress them. We use the pin/pzip system described by Buchsbaum et al.  to general optimal, contiguous partitions. For each file, we run the dynamic program on a small training set and compress the remainder of the data, the test set. Gzip results are with respect to the test sets only. Buchsbaum et al.  investigate the relationship between training size and compression performance and demonstrate a threshold after which more training data does not improve performance. Here we simply use enough training data to exceed this threshold and report this amount in Table 1. The training and test sets remain disjoint to support the validity of using a partition from a small amount of training data on a larger amount of subsequent data. In a real application, the training data would also be compressed.

All experiments were performed on one 250 MHz R10000 processor in a 24-processor SGI Challenge, with 14 GB of main memory. Each time reported is the medians of five runs.

### 4.2 Greedy Algorithms

Our hypothesis that GREEDY and GREEDYT produce partitions close in cost to that of the optimal, contiguous partition, if true implies that we can substitute the greedy algorithms for the dynamic program (DP) in purely on-line applications that cannot afford off-line training time. We thus compare compression rates of GREEDY and GREEDYT against DP and gzip, to assess the quality of the partitions; and we compare the time taken by GREEDY and GREEDYT (partitioning and compression) against gzip, to assess tractability. Table 2 shows the resulting compressed sizes using partitions computed with GREEDY and GREEDYT. Table 3 gives the time results.

GREEDY compresses to within 2% of DP on seven of the files, including four of the genetic files. It is never more than 9% bigger than DP, and with the exception of backPQQ, always outperforms gzip. GREEDYT comes within 10% of DP on seven files, including four genetic files and outperforms gzip except on backPQQ and cytoB. Both GREEDY and GREEDYT seem to outperform DP on callagen, although this would seem theoretically impossible. It is an artifact of the training/testing paradigm: we compress data distinct from that used to build the partitions.

Tables 2 and 3 show that in many cases, the greedy algorithms provide significant extra compression at acceptable time penalties. For the non-genetic files, greedy partitioning compression is less than $`1.7`$ slower than gzip yet provides 35–55% more compression. For the genetic files, the slowdown is a factor of 3–8, and the extra compression is 5–20% (ignoring backPQQ). Thus, the greedy algorithms provide a good on-line heuristic for improving compression.

### 4.3 Reordering via TSP

Our hypothesis that tour weight and compression are correlated implies that generating a TSP tour (or approximation) would yield an optimal (or near optimal) partition. Although we do not know what the optimal partition is for our files, we can assess the correlation by generating a sequence of tours and, for each, measuring the resulting compression. We also compare the compression using the best partition from the sequence against that using DP on the original ordering, to gauge the improvement yielded by reordering.

For each file, we computed various tours on the corresponding graph $`G()`$. We computed a close approximation to a TSP tour using a variation of Zhang’s branch-and-bound algorithm , discussed by Cirasella et al. . We also computed a 3-opt local optimum tour; and we used a 4-opt heuristic to compute a sequence of tours of various costs. Each tour induced an ordering of the columns. For each column ordering, we computed the optimal, contiguous partition by DP, except that we used GREEDYT on the orderings for census, due to computational limitations. Figures 1 and 2 plot the results.

The plots demonstrate a strong, positive correlation between tour cost and compression performance. In particular, each plot shows that the least-cost tour (produced by Zhang’s algorithm) produced the best compression result. Table 4 details the compression improvement from using the Zhang ordering. In five files, Zhang gives an extra compression improvement of at least 5% over DP on the original order; for cytoB, the improvement is 20%. That the original order for network outperforms the Zhang ordering is again an artifact of the training/test paradigm. Figure 1 shows that the tour-cost/compression-performance correlation remains strong for this file.

Table 5 displays the time spent computing Zhang’s tour for each file. This time is negligible compared to the time to compute the optimal, contiguous partition via DP. (The DP time on census is 168531 seconds, four orders of magnitude larger. For cytoB, the DP time is 8640 seconds, an order of magnitude larger.) Table 5 also shows that Zhang’s tour always had cost close to the Held-Karp lower bound  on the cost of the optimum TSP tour.

For off-line training, therefore, it seems that computing a good approximation to the TSP reordering before partitioning contributes significant compression improvement at minimal time cost. Furthermore, the correlation between tour cost and compression behaves similarly to what the theory in Section 2.2 would predict if $`H_𝒞()`$ were sub-additive, which suggests the existence of some other, similar structure induced by $`H_𝒞()`$ that would control this relationship.

## 5 Complexity of Table Compression: A General Framework

We now introduce a framework for studying the computational complexity of several versions of table compression problems. We start with a basic problem of finding an optimal arrangement of a set of strings to be compressed. Given a set of strings, we wish to compute an order in which to catenate the strings into a superstring $`X`$ so as to minimize the cost of compressing $`X`$ using a fixed compressor $`𝒞`$. To isolate the complexity of finding an optimal order, we restrict $`𝒞`$ to prevent it from reordering the input itself.

Let $`x=\sigma _1\mathrm{}\sigma _n`$ be a string over some alphabet $`\mathrm{\Sigma }`$, and let $`𝒞(x)`$ denote the output of $`𝒞`$ when given input $`x`$. We allow $`𝒞`$ arbitrary time and space, but we require that it process $`x`$ monotonically. That is, it reads the symbols of $`x`$ in order; after reading each symbol, it may or may not output a string. Let $`𝒞(x)_j`$ be the catenation of all the strings output, in order, by $`𝒞`$ after processing $`\sigma _1\mathrm{}\sigma _j`$. If $`𝒞`$ actually outputs a (non-null) string after reading $`\sigma _j`$, then we require that $`𝒞(x)_j`$ must be a prefix of $`𝒞(\sigma _1\mathrm{}\sigma _jy)`$ for any suffix $`y`$. We assume a special end-of-string character not in $`\mathrm{\Sigma }`$ that implicitly terminates every input to $`𝒞`$.

Intuitively, this restriction precludes $`𝒞`$ from reordering its input to improve the compression. Many compression programs used in practice work within this restriction: e.g., gzip and compress.

We use $`|𝒞(x)|`$ to abstract the length of $`𝒞(x)`$. A common measure is bits, but other measures are more appropriate in certain settings. For example, when considering LZ77 compression , we will denote by $`|𝒞(x)|`$ the number of phrases in the LZ77 parsing of $`x`$, which suffices to capture the length of $`𝒞(x)`$ while ignoring technical details concerning how phrases are encoded.

Let $`X=\{x_1,\mathrm{},x_n\}`$ be a set of strings. A batch of $`X`$ is an ordered subset of $`X`$. A schedule of $`X`$ is a partition of $`X`$ into batches. A batch $`B=(x_{i_1},\mathrm{},x_{i_s})`$ is processed by $`𝒞`$ by computing $`𝒞(B)=𝒞(x_{i_1}\mathrm{}x_{i_s})`$; i.e., by compressing the superstring formed by catenating the strings in $`B`$ in the order given. A schedule $`𝒮`$ of $`X`$ is processed by $`𝒞`$ by processing its batches, one by one, in any order. While $`𝒞(𝒮)`$ is ambiguous, $`|𝒞(𝒮)|=_{B𝒮}|𝒞(B)|`$ is well defined. Our main problem can be stated as follows.

###### Problem 5.1

Let $`X`$ be a set of strings. Find a schedule $`𝒮`$ of $`X`$ minimizing $`|𝒞(𝒮)|`$ among all schedules.

The classical shortest common superstring (SCS) problem can be phrased in terms of Problem 5.1. For two strings $`x`$ and $`y`$, let $`\text{pref}(x,y)`$ be the prefix of $`x`$ that ends at the longest suffix-prefix match of $`x`$ and $`y`$. Let $`X`$ be a set of $`n`$ strings, and let $`\pi `$ be a permutation of the integers in $`[1,n]`$. Define $`S(X,\pi )=\text{pref}(x_{\pi _1},x_{\pi _2})\text{pref}(x_{\pi _2},x_{\pi _3})\mathrm{}\text{pref}(x_{\pi _{n1}},x_{\pi _n})x_{\pi _n}`$. $`S(X,\pi )`$ is a superstring of $`X`$; $`\pi `$ corresponds to a schedule of $`X`$; and the SCS of $`X`$ is $`S(X,\pi )`$ for some $`\pi `$ . Therefore, finding the SCS is an instance of Problem 5.1, where $`𝒞()`$ is $`S()`$. Since finding the SCS is MAX-SNP hard , Problem 5.1 is MAX-SNP hard in general. Different results can hold for specific compressors, however.

We now formalize table compression problems in this framework. Consider a table $`T`$ with $`m`$ rows and $`n`$ columns, each entry a symbol in $`\mathrm{\Sigma }`$. Let $`T^c`$ be the string formed by catenating the columns of $`T`$ in order; let $`T^r`$ be the string formed by catenating the rows of $`T`$ in order.

We view $`T`$ as a set of columns $`\{T[1],\mathrm{},T[n]\}`$. A batch $`B=(T[i_1],\mathrm{},T[i_s])`$ then corresponds to a table $`T_B=T[i_1]\mathrm{}T[i_s]`$, which we can compress in column- or row-major order. A column-major order schedule $`𝒮^c`$ of $`T`$ has compression cost $`|𝒮^c|=_{B𝒮^c}|𝒞(T_B^c)|`$. A row-major order schedule $`𝒮^r`$ of $`T`$ has compression cost $`|𝒮^r|=_{B𝒮^r}|𝒞(T_B^r)|`$.

###### Problem 5.2

Given a table $`T`$, find a column-major schedule $`𝒮^c`$ of $`T`$ minimizing $`|𝒮^c(T)|`$ among all such schedules.

###### Problem 5.3

Given a table $`T`$, find a row-major schedule $`𝒮^r`$ of $`T`$ minimizing $`|𝒮^r(T)|`$ among all such schedules.

In either column- or row-major order, batches of $`T`$ are subsets of columns. In column-major order, each column of $`T`$ remains a distinct substring in any schedule. In row-major order, however, the individual strings that form a schedule are the row-major renderings of batches of $`T`$. This distinction is subtle yet crucial. Problem 5.2 becomes equivalent to Problem 5.1, so we may consider the latter in order to establish lower bounds for the former. Problem 5.3, however, is not identical to problem 5.1: the row-major order rendering of the batches results in input strings being intermixed. We emphasize this distinction in Section 7, where we show that, when $`𝒞`$ is run length encoding, Problem 5.2 can be solved in polynomial time, while Problem 5.3 is MAX-SNP hard. The connection between table compression and SCS through Problem 5.1 makes these problems theoretically elegant as well as practically motivated.

## 6 Complexity with LZ77

We use the standard definitions of L-reduction and MAX-SNP . Let $`A`$ and $`B`$ be two optimization (minimization or maximization) problems. Let $`\text{cost}_A(y)`$ be the cost of a solution $`y`$ to some instance of $`A`$; let $`\text{opt}_A(x)`$ be the cost of an optimum solution for an instance $`x`$ of $`A`$; and define analogous metrics for $`B`$. $`A`$ L-reduces to $`B`$ if there are two polynomial-time functions $`f`$ and $`g`$ and constants $`\alpha ,\beta >0`$ such that:

* Given an instance $`a`$ of $`A`$, $`f(a)`$ is an instance of $`B`$ such that $`\text{opt}_B(f(a))\alpha \text{opt}_A(a)`$;
* Given a solution $`y`$ to $`f(a)`$, $`g(y)`$ is a solution to $`a`$ such that $`|\text{cost}_A(g(y))\text{opt}_A(a)|\beta |\text{cost}_B(y)\text{opt}_B(f(a))|`$.

The composition of two L-reductions is also an L-reduction. A problem is MAX-SNP hard  if every problem in MAX-SNP can be L-reduced to it. If $`A`$ L-reduces to $`B`$, then if $`B`$ has a polynomial-time approximation scheme (PTAS), so does $`A`$. A MAX-SNP hard problem is unlikely to have a PTAS .

Now recall the LZ77 parsing rule , which is used by compressors like gzip. Consider a string $`z`$, and, if $`|z|1`$, let $`z^{}`$ denote the prefix of $`z`$ of length $`|z|1`$. If $`|z|2`$, then define $`z^{}=(z^{})^{}`$.

LZ77 parses $`z`$ into phrases, each a substring of $`z`$. Assume that LZ77 has already parsed the prefix $`z_1\mathrm{}z_{i1}`$ of $`z`$ into phrases $`z_1,\mathrm{},z_{i1}`$, and let $`z^{}`$ be the remaining suffix of $`z`$. LZ77 selects the $`i`$’th phrase $`z_i`$ as the longest prefix of $`z^{}`$ that can be obtained by adding a single character to a substring of $`(z_1\mathrm{}z_{i1}z_i)^{}`$. Therefore, $`z_i`$ has the property that $`z_i^{}`$ is a substring of $`(z_1z_2\mathrm{}z_{i1}z_i)^{}`$, but $`z_i`$ is not a substring of $`(z_1z_2\mathrm{}z_{i1}z_i)^{}`$. This recursive definition is sound .

After parsing $`z_i`$, LZ77 outputs an encoding of the triplet $`(p_i,\mathrm{}_i,\alpha _i)`$, where $`p_i`$ is the starting position of $`z_i^{}`$ in $`z_1z_2\mathrm{}z_{i1}`$; $`\mathrm{}_i=|z_i|1`$; and $`\alpha _i`$ is the last character of $`z_i`$. The length of the encoding is linear in the number of phases, so when $`𝒞`$ is LZ77, we denote by $`|𝒞(z)|`$ the number of phrases in the parsing of $`z`$. This cost function is commonly used to establish the performance of LZ77 parsing .

### 6.1 Problem 5.1

We show that Problem 5.1 is MAX-SNP hard when $`𝒞`$ is LZ77. Consider TSP(1,2), the traveling salesman problem on a complete graph where each distance is either 1 or 2. An instance of TSP(1,2) can be specified by a graph $`H`$, where the edges of $`H`$ connect those pairs of vertices with distance 1. The problem remains MAX-SNP hard if we further restrict the problem so that the degree of each vertex in $`H`$ is bounded by some arbitrary but fixed constant . This result holds for both symmetric and asymmetric TSP(1,2); i.e., for both undirected and directed graphs $`H`$. We assume that $`H`$ is directed. The following lemma shows that we may also assume that no vertex in $`H`$ has outdegree 1.

###### Lemma 6.1

TSP(1,2) L-reduces to TSP(1,2) with the additional stipulation that no vertex has only one outgoing cost-1 edge.

Proof.  Consider instance $`A`$ of TSP(1,2). For each vertex $`v`$ with only one outgoing cost-1 edge, to some $`v^{}`$, we create a new vertex $`v^{\prime \prime }`$ such that edges $`(v,v^{\prime \prime })`$, $`(v^{\prime \prime },v)`$, and $`(v^{\prime \prime },v^{})`$ have cost 1 and all other edges incident on $`v^{\prime \prime }`$ have cost 2. Thus we form instance $`B`$. A solution to $`B`$ is mapped to a solution to $`A`$ by splicing out all newly created vertices. If $`A`$ has $`n`$ vertices, then $`B`$ has at most $`2n`$ vertices. All solutions to both have cost $`O(n)`$, so we need only prove that the reverse mapping of solutions preserves optimality.

Assume $`S_B`$ is an optimal solution to $`B`$ and $`S_A`$ the corresponding mapped solution to $`A`$. Note that $`\text{cost}(S_A)\text{cost}(S_B)\eta `$, where $`\eta `$ is the number of vertices created to form $`B`$. (We drop the subscripts to $`\text{cost}()`$, as there is no ambiguity.) If $`S_A`$ is not optimal, there is some $`S_A^{}`$ such that $`\text{cost}(S_A^{})<\text{cost}(S_A)`$. We can form a solution $`S_B^{}`$ to $`B`$ by replacing each edge $`(v,z)`$ in $`S_A^{}`$, where $`v`$ has only one cost-1 outgoing edge, with edges $`(v,v^{})`$ and $`(v^{},z)`$. This gives $`\text{cost}(S_B^{})=\text{cost}(S_A^{})+\eta <\text{cost}(S_A)+\eta \text{cost}(S_B)`$, contradicting the optimality of $`S_B`$. $`\mathrm{}`$

We associate a set $`S(H)`$ of strings to the vertices and edges of $`H`$; $`S(H)`$ will be the input to Problem 5.1. Each vertex $`v`$ engenders three symbols: $`v`$, $`v^{}`$, and $`\$_v`$. Let $`w_0,\mathrm{},w_{d1}`$ be the vertices on the edges out of $`v`$ in $`H`$, in some arbitrary but fixed cyclic order. For $`0i<d`$ and mod-$`d`$ arithmetic, we say that edge $`(v,w_i)`$ cyclicly precedes edge $`(v,w_{i+1})`$. The $`d+1`$ strings we associate to $`v`$ and these edges are: $`e(v,w_i)=(v^{}w_{i1})^4v^{}w_i`$, for $`0i<d`$ and mod-$`d`$ arithmetic; and $`s(v)=v^4(v^{})^5\$_v`$. That $`d1`$ implies that $`w_iw_{i+1}`$ when $`d0`$,

To prove MAX-SNP hardness, we first show how to transform a TSP(1,2) solution for $`H`$ into a solution to Problem 5.1 with input $`S(H)`$. We then show how to transform in polynomial time a solution to Problem 5.1 into a TSP(1,2) solution of a certain cost. We use the intermediate step of transforming the first solution into a canonical form of at most the same cost.

The canonical form solution will correspond to the required TSP(1,2) tour. We will show that, for all edges $`(v,w)`$, $`e(v,w)`$ will parse into one phrase when immediately preceded by $`e(v,y)`$ for the edge $`(v,y)`$ that cyclicly precedes $`(v,w)`$, and into more than one phrase otherwise; and we will show that $`s(v)`$ will parse into two phrases when immediately preceded by $`e(x,v)`$ for some edge $`(x,v)`$, and into three phrases otherwise. Thus, an edge $`(v,w_i)`$ in the path will best be encoded as $`s(v)e(v,w_{i+1})e(v,w_{i+2})\mathrm{}e(v,w_i)s(w_i)`$. This is the core idea of our canonical form.

Lemmas 6.26.6 provide a few needed facts about the parsing of strings in $`S(H)`$. In what follows, $`X`$ denotes both a batch in $`S(H)`$ and the string obtained by catenating the strings in the batch in order.

###### Lemma 6.2

Let $`X=x_1\mathrm{}x_s`$ be a batch of $`S(H)`$, where each $`x_i`$ is $`s(v)`$ for some vertex $`v`$ or $`e(v,w)`$ for some edge $`(v,w)`$. For each $`1js`$, some phrase in the LZ77 parsing of $`X`$ ends at the last symbol of $`x_j`$.

Proof.  The proof is by induction. The base case is for $`j=1`$. If $`x_1=s(v)`$ for some vertex $`v`$, then the lemma holds, because $`\$_v`$ appears only at the end of $`s(v)`$. Otherwise, $`x_1=e(v,w_i)`$ for some edge $`(v,w_i)`$. Since $`x_1`$ appears first in $`X`$, its parsing is $`v^{}`$, $`w_{i1}`$, $`(v^{}w_{i1})^3v^{}w_i`$. (That no vertex has outdegree one implies that $`w_iw_{i1}`$.) The lemma again holds.

Assume by induction that the lemma is true up through the parsing of $`x_{j1}`$; we show that it holds for the parsing through $`x_j`$. Again, if $`x_j=s(v)`$ for some vertex $`v`$, the lemma is true, because $`\$_v`$ appears only at the end of $`s(v)`$. Otherwise, $`x_j=e(v,w_i)`$, for some edge $`(v,w_i)`$. There are two cases.

1. $`x_{j1}=e(v,w_{i1})`$. Then $`x_{j1}x_j=(v^{}w_{i2})^4(v^{}w_{i1})^5v^{}w_i`$. By induction, a phrase ends at the first occurrence of $`w_{i1}`$. Thus, the next phrase is $`(v^{}w_{i1})^4v^{}w_i=x_j`$.
2. $`x_{j1}e(v,w_{i1})`$. Again by induction, the first phrase, say $`c`$, of the parsing that overlaps $`x_j`$ must start at the first character of $`x_j`$. Since $`(v^{}w_{i1})^2`$ does not occur in $`x_1\mathrm{}x_{j1}`$, the first phrase cannot extend past the fourth character of $`x_j`$. We have the following subcases.

1. $`c`$ ends at the first character of $`x_j`$. Therefore $`v^{}`$ does not occur in $`x_1\mathrm{}x_{j1}`$. Since $`x_j=(v^{}w_{i1})^4v^{}w_i`$, we have that the phrase following $`c`$, say $`c^{}`$, must be either $`w_{i1}`$ or $`w_{i1}v^{}`$, depending on whether or not $`w_{i1}`$ occurs in $`x_1\mathrm{}x_{j1}`$. (1) When $`c^{}`$ is $`w_{i1}`$, the next phrase is $`(v^{}w_{i1})^3v^{}w_i`$ and ends on the last character of $`x_j`$, as required. (2) When $`c^{}`$ is $`w_{i1}v^{}`$, the next phrase is $`(w_{i1}v^{})^3w_i`$, again completing the induction.
    2. Remaining cases. When $`c`$ ends at the second and third character of $`x_j`$, the result follows as in (2a.1) and (2a.2), respectively. When $`c`$ ends at the forth character, the next phrase is $`(v^{}w_{i1})^2v^{}w_i`$ and ends at the last character of $`x_j`$ as required.

$`\mathrm{}`$

###### Lemma 6.3

Let $`X`$ be a batch of $`S(H)`$ and $`v`$ be any vertex such that $`s(v)X`$. If $`s(v)`$ is immediately preceded by $`e(q,v)`$ for some edge $`(q,v)`$, $`s(v)`$ is parsed into precisely two phrases during the parsing of $`X`$; otherwise, $`s(v)`$ is parsed into precisely three phrases.

Proof.  Assume first that $`s(v)`$ is immediately preceded by $`e(q,v)`$ for some edge $`(q,v)`$. Then $`e(q,v)s(v)=(q^{}z)^4q^{}vv^4(v^{})^5\$_v`$ for some $`z`$. By Lemma 6.2, a phrase of the parsing must end with the last character of $`e(q,v)`$. Since $`v^4`$ does not appear elsewhere in $`X`$, the next two phrases of the parsing must be $`v^4v^{}`$ and $`(v^{})^4\$_v`$.

In the other case, $`v^2`$ does not occur to the left of $`s(v)`$ in $`X`$. Again using Lemma 6.2, the parsing of $`X`$ has a phrase starting at $`s(v)`$. If $`v`$ appears to the left of $`s(v)`$ in $`X`$, the parsing produces $`v^2`$, $`v^2v^{}`$, and $`(v^{})^4\$_v`$; Otherwise, it produces $`v`$, $`v^3v^{}`$, and $`(v^{})^4\$_v`$. $`\mathrm{}`$

###### Lemma 6.4

Let $`X`$ be a batch of $`S(H)`$ and $`(v,w)`$ be any edge such that $`e(v,w)X`$. Let $`(v,y)`$ be the edge that cyclicly precedes $`(v,w)`$. If $`e(v,w)`$ is immediately preceded in $`X`$ by $`e(v,y)`$, then $`e(v,w)`$ is parsed into precisely one phrase during the parsing of $`X`$; if $`e(v,w)`$ is immediately preceded by $`s(v)`$, then $`e(v,w)`$ is parsed into precisely two phrases; in any other case, $`e(v,w)`$ is parsed into at least two phrases.

Proof.  By Lemma 6.2, some phrase starts at the first character of $`e(v,w)`$. Assume $`e(v,y)`$ immediately precedes $`e(v,w)`$; $`e(v,y)e(v,w)=(v^{}z)^4v^{}y(v^{}y)^4v^{}w`$ for some $`z`$. The parsing of $`e(v,w)`$ produces the one phrase $`(v^{}y)^4v^{}w=e(v,w)`$. (Nowhere else does this string appear in $`X`$.)

Assume $`s(v)`$ immediately precedes $`e(v,w)`$; $`s(v)e(v,w)=v^4(v^{})^5\$_v(v^{}y)^4v^{}w`$. If $`v^{}y`$ occurs earlier in $`X`$, the parsing of $`e(v,w)`$ produces phrases $`v^{}yv^{}`$ and $`(yv^{})^3w`$, because $`v^{}yv^{}`$ cannot occur elsewhere. Otherwise, the parsing produces $`v^{}y`$ and $`(v^{}y)^3v^{}w`$.

In any other case, $`e(v,w)`$ is preceded by a character other than $`v^{}`$. If $`v^{}`$ occurs earlier in $`X`$, then the parsing of $`e(v,w)`$ produces two phrases as in the case of $`s(v)`$ preceding $`e(v,w)`$. Otherwise, the parsing produces $`v^{}`$ and then at least one more phrase. $`\mathrm{}`$

Now define a schedule $`Y_1,\mathrm{},Y_t`$ to be standard if and only if: for each batch $`Y_i`$, the order in which the strings $`s(v)`$ appear in $`Y_i`$ corresponds to a path in $`H`$; the paths associated to $`Y_i`$ and $`Y_j`$ are disjoint for each $`ij`$; and each vertex of $`H`$ appears as $`s(v)`$ in some batch $`Y_i`$.

We give a polynomial time algorithm that transforms a schedule $`𝒮=(X_1,\mathrm{},X_g)`$ into a standard schedule that parses into no more phrases than does $`𝒮`$. The algorithm consists of two phases. The first phase computes a set of disjoint paths that covers all the vertices of $`H`$. It iteratively combines paths, guided by $`𝒮`$, until no further combination is possible. The second phase transforms each path into a batch such that the resulting schedule is standard.

Algorithm STANDARD

* 1. Place each vertex $`v`$ of $`H`$ in a single-vertex path. If $`s(v)`$ is the first string in some batch in $`𝒮`$, label $`v`$ terminal; otherwise, label $`v`$ nonterminal.
    2. While there exists a path with nonterminal left end point, pick one such end point $`v`$ and process it as follows. Let $`X_i`$ be the batch in which $`s(v)`$ occurs. Let $`x(u)`$ be the string (associated to either vertex $`u`$ or to one of its outgoing edges) that precedes $`s(v)`$ in $`X_i`$. If $`x(u)`$ ends in a symbol other than $`v`$, label vertex $`v`$ terminal. Otherwise, $`(u,v)H`$, so connect $`u`$ to $`v`$, and, for each edge $`(u,w)H`$, $`uw`$, such that $`s(w)`$ is immediately preceded by $`e(u,w)`$, declare $`w`$ terminal. (This guarantees that Phase One actually builds paths.)
* Let $`𝒜_1,\mathrm{},𝒜_t`$ be the paths obtained at the end of Phase One. We transform each path $`𝒜_j`$ into a batch $`Y_j`$. If $`𝒜_j`$ consists of a single vertex $`v`$, then $`Y_j`$ consists of $`s(v)`$ followed by all the $`e(v,w_j)`$’s arranged in cyclic order.

Otherwise, $`𝒜_j`$ contains more than one vertex. Initially $`Y_j`$ is empty. For each edge $`(u,v)`$ in order in the path, we append to $`Y_j`$: $`s(u)`$ followed by all of its $`e(u,w_j)`$’s, in cyclic order ending with $`e(u,v)`$. When there are no more edges to process, the last vertex of the path is processed as in the singleton-vertex case.

###### Lemma 6.5

In polynomial time, Algorithm STANDARD transforms schedule $`X_1,\mathrm{},X_g`$ into a standard schedule $`Y_1,\mathrm{},Y_t`$ of no higher cost.

Proof.  That Algorithm STANDARD runs in polynomial time and $`Y_1,\mathrm{},Y_t`$ is standard follow immediately from the specification.

We now show that each batch $`Y_j`$ parses into no more phrases than do its corresponding components in the input schedule. Consider the path, $`(v_1,v_2,\mathrm{},v_r)`$ from which $`Y_j`$ is derived. Let $`d(v)`$ be the outdegree of any vertex $`v`$. $`Y_j=s(v_1)e(v_1,w_1^1)\mathrm{}e(v_1,w_{d(v_1)}^1)\mathrm{}s(v_r)e(v_r,w_1^r)\mathrm{}e(v_r,w_{d(v_r)}^r)`$, where the $`w_j^i`$’s are the neighbors in cyclic order out of $`v_i`$ and, for $`1i<r`$, we assume without loss of generality that $`w_{d(v_i)}^i=v_{i+1}`$.

By Lemma 6.3, for $`2ir`$, $`s(v_i)`$ parses into two phrases, which is optimal. By Lemma 6.4, for $`1ir`$ and $`2jd(v_i)`$, $`e(v_i,w_j^i)`$ parses into one phrase, which is optimal. We thus need only consider the parsing of $`s(v_1)`$ and, for $`1ir`$, $`e(v_i,w_1^i)`$.

The strings $`e(v_i,w_1^i)`$, for $`1ir`$ each parse into two phrases in $`Y_j`$, by Lemma 6.4. There must be some $`e(v_i,x)`$ that is not immediately preceded by its cyclic predecessor in some $`X_k`$, and this instance of $`e(v_i,x)`$ also parses into at least two phrases, by Lemma 6.4. This accounts for the first $`e()`$ string immediately following each $`s()`$ string in $`Y_j`$.

Finally, if $`s(v_1)`$ is not immediately preceded by some $`e(v,v_1)`$ in the input batch $`X_k`$ in which $`s(v_1)`$ appears, we are done, for $`s(v_1)`$ is parsed into three phrases in both $`X_k`$ and $`Y_j`$, by Lemma 6.3. Otherwise, consider the maximal sequence $`e(v,w_a)e(v,w_{a+1})\mathrm{}e(v,w_{a+\mathrm{}}=v_1)s(v_1)`$ in $`X_k`$, where the $`w_a`$’s are cyclicly ordered neighbors of $`v`$. Because STANDARD declared $`v_1`$ to be terminal, there was another edge $`(v,y)`$ such that $`e(v,y)`$ immediately preceded $`s(y)`$ in some $`X_k^{}`$, which STANDARD used to connect $`v`$ and $`y`$ in some path. This engenders an analogous maximal chain of $`e(v,)`$ strings followed by $`s(y)`$ in $`X_k^{}`$.

Thus, there are at least two strings $`e(v,)`$ not immediately preceded in the input by their cyclic predecessors; Lemma 6.4 implies each is parsed into at least two phrases. We can charge the extra phrase generated by $`s(v_1)`$ in $`Y_j`$ against one of them, leaving the other for the extra phrase in the parsing of the $`e(v,)`$ phrase immediately following $`s(v)`$ in some $`Y_j^{}`$. $`\mathrm{}`$

###### Lemma 6.6

A batch $`Y_j`$ output by STANDARD, corresponding to a path $`(v_1,\mathrm{},v_r)`$, parses into $`3r+1+_{i=1}^rd(v_i)`$ phrases.

Proof.  By Lemma 6.3, each $`s()`$ string parses into 2 phrases, except $`s(v_1)`$, which parses into 3, contributing $`2r+1`$ phrases. Lemma 6.4 implies that each $`e()`$ parses into 1 phrase, except each following an $`s()`$, which parses into 2, contributing $`r+_{i=1}^rd(v_i)`$ phrases. $`\mathrm{}`$

###### Theorem 6.7

Problem 5.1 is MAX-SNP hard when $`𝒞`$ is LZ77.

Proof.  Let the graph $`H`$ defined at the beginning of the section have $`n_h`$ vertices and $`m_h`$ edges. Let $`k`$ be the minimum number of cost-2 edges that suffice to form a TSP(1,2) solution. Then the cost of the solution is $`n_h1+k`$. Associating strings to vertices and edges of $`H`$, as discussed above, we argue that the optimal schedule for those strings produces $`m_h+k+3n_h+1`$ phrases. The reduction is linear, since $`m_h=O(n_h)`$ by the assumption of bounded outdegree.

Assume that the TSP(1,2) solution with $`k`$ cost-2 edges is the path $`v_1,v_2,\mathrm{},v_{n_h}`$. Then in polynomial time we can construct a corresponding standard schedule of the form output by STANDARD, which Lemma 6.6 shows parses into $`m_h+k+3n_h+1`$ phrases.

For the converse, assume that we are given a schedule of cost $`m_h+k+3n_h+1`$. By Lemma 6.5, we can transform it in polynomial time into a standard schedule $`Y_1,\mathrm{},Y_t`$ of no higher cost. Recall that to each batch we can associate a path of $`H`$. Let $`v_1,v_2,\mathrm{},v_{n_h}`$ the an ordering of the vertices of $`H`$ corresponding to an arbitrarily chosen processing order for the sequence of batches. Then, $`H`$ cannot be missing more than $`k`$ of the edges $`(v_i,v_{i+1})`$, or else, by Lemma 6.6, the cost of the standard schedule would exceed $`m_h+k+3n_h+1`$. $`\mathrm{}`$

## 7 Complexity with Run Length Encoding

In run length encoding (RLE), an input string is parsed into phrases of the form $`(\sigma ,n)`$, where $`\sigma `$ is a character, and $`n`$ is the number of times $`\sigma `$ appears consecutively. For example, $`aaaabbbbaaaa`$ is parsed into $`(a,4)(b,4)(a,4)`$.

### 7.1 Problem 5.1

###### Theorem 7.1

Problem 5.1 can be solved in polynomial time when $`𝒞`$ is run length encoding.

Proof.  Let $`x_1,\mathrm{},x_n`$ be the input strings. We can assume without loss of generality that each $`x_i`$ is of the form $`\sigma \sigma ^{}`$; i.e., two distinct characters. The parsing of any characters between them cannot be optimized by rearranging the strings. Furthermore, if $`x_i=\sigma \sigma `$, we can simply merge $`x_i`$ with another string, $`x_j`$, that begins or ends with $`\sigma `$; if no such $`x_j`$ exists, we can ignore $`x_i`$ completely, since again its parsing cannot be optimized by rearrangement.

We claim that a shortest common superstring (SCS) of the input corresponds to an optimal schedule. As described earlier, an SCS is $`\text{pref}(x_{\pi _1},x_{\pi _2})\mathrm{}\text{pref}(x_{\pi _{n1}},x_{\pi _n})x_{\pi _n}`$ for some permutation $`\pi `$. Note that $`\text{pref}(x_i,x_j)`$ is of length 2 if the last character of $`x_i`$ equals the first of $`x_j`$ and 3 otherwise. Thus, an SCS gives an optimal RLE parsing, and SCS can be solved in polynomial time when all input strings are of length two . $`\mathrm{}`$

### 7.2 Problem 5.3

As in Section 6.1, we transform the vertices and edges of $`H`$ into an instance of Problem 5.3. We associate a column to each vertex and edge of $`H`$.

For each vertex $`v`$, we generate three symbols: $`v`$, $`v^{}`$, and $`v^{\prime \prime }`$. Let $`w_0,\mathrm{},w_{d1}`$ be the vertices on the edges out of $`v`$ in some fixed, arbitrary cyclic order. We associate the following strings to $`v`$ and its outgoing edges: $`s(v)=v^{}v^{\prime \prime }v`$; and $`e(v,w_i)=v^{}v^{\prime \prime }w_i`$, $`0i<d`$. The input table is formed by assigning each such string, over all the vertices, to a column.

Consider a TSP(1,2) solution with $`k`$ cost-2 edges. We can arrange the induced strings into a table $`T`$ describing these paths. Place all strings corresponding to a vertex $`v`$ in a contiguous interval of the table with $`s(v)`$ being the first column of the interval. For any edge $`(v,q)`$ in the collection of paths, place the interval corresponding to $`q`$ immediately after that corresponding to $`v`$, and place the string $`e(v,q)`$ last in the interval for $`v`$; otherwise, the order of the intervals and of the strings corresponding to edges can be arbitrary. We say the table is in standard form for the collection of paths.

###### Theorem 7.2

Problem 5.3 is MAX-SNP hard for tables of at least 3 rows when $`𝒞`$ is run length encoding.

Proof.  We prove the theorem for three rows first and then extend it to larger numbers of rows. Let $`n_h`$ and $`m_h`$ be the number of vertices and edges in $`H`$, rsp., and let $`n=n_h+m_h`$ be the number of columns in the induced table. Associate strings to the vertices and edges as described above. Let $`k`$ be the minimum number cost-2 edges that suffice to form a TSP(1,2) solution for $`H`$. Then the cost of the solution is $`n_h1+k`$. Let $`v_1,v_2,\mathrm{},v_{n_h}`$ be an ordering of the vertices in $`H`$ corresponding to the $`k+1`$ disjoint paths. Let $`T`$ be the corresponding standard form table. Let $`S`$ be the schedule obtained by taking as a batch each interval of the table corresponding to a path. The row-major cost of $`S`$ is $`2n+m_h+k+1`$. This completes one direction of the transformation.

As for the other direction, assume that we are given a solution to the instance of optimum table compression that has cost $`2n+m_h+k+1`$. Let $`T^{}`$ be the table of the solution schedule. In polynomial time, we can transform $`T^{}`$ into a standard form table $`T`$ with a schedule of at most the same cost. We simply observe that, if the $`e()`$ and $`s()`$ strings for any vertex are not contiguous, we can rearrange the columns to make them so, saving at least two phrases and generating at most two in the new parsing.

Since a table in standard form corresponds to an ordering of the vertices, it must be that $`H`$ cannot be missing more than $`k`$ edges, or else the cost of the table in standard form would be greater than $`2n+m_h+k+1`$.

When the number of rows $`m`$ exceeds three, we use one additional character $`\$`$. Each string is as in the case $`m=3`$, except that now is augmented to end with the suffix $`\$^{m3}`$. This would add one more phrase to the parsing of the set of strings, and the linearity of the transformation still holds. $`\mathrm{}`$

## 8 Conclusion

We demonstrate a general framework that links independence among groups of variables to efficient partitioning algorithms. We provide general solutions in ideal cases in which dependencies form equivalence classes or cost functions are sub-additive. The application to table compression suggests that there also exist weaker structures that allow partitioning to produce significant cost improvements. Open is the problem of refining the theory to explain these structures and extending it to other applications.

Based on experimental results, we conjecture that our TSP reordering algorithm is close to optimal; i.e., that no partition-based algorithm will produce significantly better compression rates. It is open if there exists a measurable lower bound for compression optimality, analogous, e.g., to the Held-Karp TSP lower bound.

Finally, while we have shown some MAX-SNP hardness results pertaining to table compression, it is open whether the problem is even approximable to within constant factors.

## Acknowledgements

We are indebted to David Johnson for running his implementation of Zhang’s algorithm and local 3-opt on our files. We thank David Applegate, Flip Korn, Cecilia LaNave, S. Muthukrishnan, Grazieno Pesole, and Andrea Sgarro for many useful discussions.


format: text

text:
# Model of coarsening and vortex formation in vibrated granular rods

\[

## Abstract

Neicu and Kudrolli observed experimentally spontaneous formation of the long-range orientational order and large-scale vortices in a system of vibrated macroscopic rods. We propose a phenomenological theory of this phenomenon, based on a coupled system of equations for local rods density and tilt. The density evolution is described by modified Cahn-Hilliard equation, while the tilt is described by the Ginzburg-Landau type equation. Our analysis shows that, in accordance to the Cahn-Hilliard dynamics, the islands of the ordered phase appear spontaneously and grow due to coarsening. The generic vortex solutions of the Ginzburg-Landau equation for the tilt correspond to the vortical motion of the rods around the cores which are located near the centers of the islands.

\]

Vibrated granular materials exhibit many interesting phenomena, including formation of cellular and localized patterns, convection, phase separation etc. . Recently, Neicu and Kudrolli  studied the dynamics of a layer of long cylindrical grains (rods) subjected to vertical vibration, and discovered a surprising phenomenon of spontaneous formation of islands of vertically aligned rods which co-exit with the “sea” of randomly packed almost horizontal rods. Subsequently, small islands merge and form large islands which typically exhibit collective vortical motion of rods. Near the core of the vortex this motion has a form of solid body rotation, while farther away from the core, the angular velocity decays.

While the bistability and first order phase transitions leading to phase separation and coarsening are typical in granular dynamics (observed, for example, in monolayer excitation, vibrated powder heaping , electrostatically-driven granular systems ) and usually caused by the inelasticity of grains, the emergence of the vortical motion within the ordered phase is rather unexpected. Experiment  shows that rods within a vortex are tilted in the azimuthal direction, and slowly drift in the direction of the tilt. Authors  suggested that the drift occurs due to the confinement of the rod vibration by its tilted neighbors.

In this Letter we introduce a continuous phenomenological model of the transition to the ordered vortical state based on the modified Cahn-Hilliard equation governing the dynamics of local rods density and the Ginzburg-Landau type equation for the tilt. Our model reproduces qualitatively the observed phase separation, coarsening and vortex formation. We derive the solutions for the stationary vortices and discuss their stability.

Model. The motion of rods is described by the momentum conservation equation in the form

$$\rho \left(\frac{D𝐯}{Dt}+\zeta 𝐯\right)=p+\alpha 𝐧f_0(n)\rho $$
(1)
Here $`𝐯=(v_x,v_y)`$ is the horizontal velocity of rods, $`\rho `$ is the density, $`p`$ is the hydrodynamic pressure, the tilt vector $`𝐧=(n_x,n_y)`$ is the projection of the rod on $`xy`$ plane, and $`n=|𝐧|`$. The term $`\alpha 𝐧f_0(n)\rho `$, accounts for the driving force from the vibrating bottom exerted on the tilted rod. According to experiments , the driving force is proportional to the tilt of rods for small tilt values, but the saturates and eventually decays to zero at $`n>n_0`$. Finally, the term $`\zeta 𝐯`$ describes the momentum dissipation due to bottom friction. Eq. (1) must be augmented by the mass conservation equation in the from

$$_t\rho +\mathrm{div}(𝐯\rho )=0$$
(2)
In the following we assume that the friction is strong, and neglect the inertia term $`\frac{D𝐯}{Dt}`$ with respect to the friction term $`\zeta 𝐯`$. Thus, we can express the velocity in the form

$$𝐯=\frac{1}{\zeta \rho }\left(p\alpha 𝐧f_0(n)\rho \right)$$
(3)
Now, substituting the velocity (3) into the mass conservation law Eq. (2), we obtain

$$_t\rho =\zeta ^1\mathrm{div}\left(p\alpha 𝐧f_0(n)\rho \right)$$
(4)
To describe phenomenologically the experimentally observed phase separation and coarsening we employ the Cahn-Hilliard approach (see for review ). We assume that the pressure $`p`$ can be obtained from the variation of certain “free energy” type functional of the $`\rho `$ field

$$p=\frac{\delta F}{\delta \rho }$$
(5)
We adopt the standard form of the free energy taking into account the local dynamics and diffusive-type coupling

$$F=𝑑x𝑑y\left(l^2(\rho )^2+f(\rho )\right)$$
(6)
where $`l`$ is the characteristic length scale related to the rod thickness. To account for the bistability and phase separation, function $`f`$ should have two minima separated by a maximum. For simplicity we choose a quartic polynomial form of $`f`$, and define $`df/d\rho =A(\rho \rho _0)(\rho _{}\rho )(\delta \rho )`$, where $`0\rho _0<\delta <\rho _{}`$ are characteristic rod densities depending on the driving acceleration . Substituting (5) and (6) into Eq. (4), after appropriate rescaling we obtain the modified Cahn-Hilliard equation

$`_t\rho `$ $`=`$ $`^2\left(^2\rho \rho (1\rho )(\delta \rho )\right)`$ (7)
$``$ $`\alpha \mathrm{div}\left(𝐧f_0(n)(\rho +\rho _0)\right),`$ (8)
where we keep the same notations for the rescaled variables and parameters.

To close the description we need to add an equation for the evolution of tilt $`𝐧`$. For density $`\rho <1`$ the vertical orientation of rods corresponding to $`𝐧=0`$ is unstable, as rods spontaneously tilt. We assume that the growth rate of the instability depends on the rods packing density, so we can write for the local dynamics

$$_t𝐧=f_1(\rho )𝐧|𝐧|^2𝐧$$
(9)
with $`f_1(\rho )=a_0a_1\rho `$, $`a_{0,1}>0`$ some constants. As $`\rho `$ increases, the instability saturates and the equilibrium tilt diminishes. In addition, rods interact with each other, which leads to the additional spatial derivative operator $`\widehat{D}[𝐧]`$ in Eq.(9). Since the tilt field is not divergence-free, from the general symmetry considerations, in the lowest (second) order, the “diffusion” operator acting on $`𝐧`$, takes the form $`\widehat{D}[𝐧]=f_2(\rho )\left(\xi _1^2𝐧+\xi _2\mathrm{div}𝐧\right)`$. The coefficients $`\xi _{1,2}`$ in this expression are analogous to the first and second viscosity in ordinary fluids (see ). Function $`f_2(\rho )`$ describes the decrease of the spatial coupling strength as the rods density decreases. We assume that in the gas phase ($`\rho 0`$) the spatial coupling between the rods is small and their tilt becomes large and uncorrelated. Accordingly, we set $`f_2=\rho `$, if $`\rho >0`$ and $`f_2=0`$ otherwise. Finally, we include the simplest term describing coupling between the tilt and the density gradient $`\beta \rho `$. Combining all these terms we arrive at

$`_t𝐧`$ $`=`$ $`f_1(\rho )𝐧|𝐧|^2𝐧+`$ (10)
$`+`$ $`f_2(\rho )\left(\xi _1^2𝐧+\xi _2\mathrm{div}𝐧\right)+\beta \rho .`$ (11)
It is convenient to introduce new complex variable $`\psi =n_x+in_y`$. Then, Eq. (11) assumes the form of the generalized Ginzburg-Landau equation ($`\overline{\xi }=\xi _1+\xi _2/2`$)

$`_t\psi `$ $`=`$ $`\left(f_1(\rho )|\psi |^2\right)\psi +\beta (_x+i_y)\rho `$ (12)
$`+`$ $`f_2(\rho )\left(\overline{\xi }^2\psi +{\displaystyle \frac{\xi _2}{2}}(_x+i_y)^2\psi ^{}\right)`$ (13)
Phase separation. Eq. (8) exhibits phase separation only in a certain range of initial conditions $`\rho `$ and parameter $`\delta `$. The total mass conservation yields $`\rho 𝑑x𝑑y=S\mathrm{\Phi }`$, where $`S`$ is the cavity area, and $`\mathrm{\Phi }`$ is the average density determined by the filling fraction. Stationary solution to Eq. (8) obeys

$$^2\rho \rho (1\rho )(\delta \rho )=B,$$
(14)
where $`B=const`$ is determined below. In the phase separation regime spatially-homogeneous Eq. (14) has three real roots $`\rho _1<\rho _2<\rho _3`$. The root $`\rho _2`$ corresponds to the unstable solution and $`\rho _{1,3}`$ to the stable ones. Final stage of the phase separation leads to one domain of the high-density phase (solid) $`\rho =\rho _3`$ of the area $`S_h`$ and the low-density one (gas) $`\rho =\rho _1`$ of the area $`S_l=SS_h`$. From the mass conservation one obtains

$$S_h\rho _3+(SS_h)\rho _1=S\mathrm{\Phi }.$$
(15)
Here we neglect the interfacial contributions assuming that $`S_l,S_h1`$. In addition, Eq. (15) must be augmented by the condition that the free energy densities of the both phases are equal, which is expressed by the relation (so-called area rule, see e.g. )

$$_{\rho _1}^{\rho _3}[\rho (1\rho )(\delta \rho )B)]d\rho =0.$$
(16)
Eq. (16) and Eq. (14) fix the value of $`B`$, and, correspondingly the roots $`\rho _{1,3}`$. ¿From Eq. (15) it follows that the solution with two phases is possible if $`\rho _1<\mathrm{\Phi }<\rho _3`$. It defines the phase separation region in Fig. 1.

Let us consider radially symmetric vortex solutions to Eq. (13). In this case $`\rho `$ is a function of the polar radius $`r`$, and $`\psi `$ can be expressed as $`\psi =\mathrm{exp}(\pm i\theta )w(r)`$, where $`w`$ is a complex function, and $`\theta `$ is the polar angle (for definiteness we take sign $`+`$). Using $`_x+i_y=\mathrm{exp}(i\theta )(_r+i/r_\theta )`$, we obtain from Eq. (13)

$$_tw=f_2(\overline{\xi }_r^2w+\frac{\xi _2}{2}_r^2w^{})+f_1w|w|^2w+\beta \rho _r$$
(17)
where $`_r^2=_r^2+r^1_rr^2`$ is the radial Laplacian operator. For $`\xi _2,\beta =0`$, Eq. (17) possesses a stationary solution in the form $`w=W(r)\mathrm{exp}(i\varphi _0)`$ with the real positive magnitude $`W`$ and an arbitrary constant phase $`\varphi _0`$. The terms $`\xi _2,\beta `$ still permit the constant phase solutions, but they destroy the continuous phase degeneracy. Indeed, Eq. (17) for $`\beta =0`$ yields

$`f_2(\overline{\xi }+{\displaystyle \frac{\xi _2}{2}}\mathrm{cos}2\varphi _0)_r^2W+f_1WW^3=0`$ (18)
and $`\mathrm{sin}2\varphi _0=0`$. Solutions exist only for $`\varphi _0=0,\pi `$ or $`\pm \pi /2`$. Function $`W`$ describes the standard (and well-documented) vortex solution to the Ginzburg-Landau equation with the property $`Wf_1^{1/2}`$ for $`r\mathrm{}`$ and $`Wr`$ for $`r1`$ (for its rational approximation see e.g. ). Solutions with $`\varphi _0=0,\pi `$ describe sinks (sources) with zero circulations, whereas the solutions with $`\varphi _0=\pm \pi /2`$ are vortices with non-zero circulations. The sign of $`\varphi _0`$ determines the direction of rotation. Near the vortex core $`Wr`$, which corresponds to the solid body rotation, as velocity $`v_\theta W`$. Far away from the core the vortex exhibits differential rotation.

For $`\beta _r\rho 0`$, Eq. (17) has constant phase solutions only with $`\varphi _0=0,\pi `$, i.e. with zero circulation. However, it does not guarantee the selection of this solution in the bulk of large islands where the density gradient is small. It is easy to show that solutions with $`\varphi _0=0`$ are energetically unfavorable with respect to rotating solutions with $`\varphi _0=\pm \pi /2`$ if $`\beta _r\rho `$ is small. Indeed, Eq. (17) for $`\beta _r\rho =0`$ can be written in the form

$$_tw=\frac{\delta U}{\delta w^{}}$$
(19)
with the free energy functional

$`U`$ $`=`$ $`{\displaystyle }d𝐫[f_2\overline{\xi }(|_rw|^2+r^2|w|^2)`$ (20)
$`+`$ $`{\displaystyle \frac{f_2\xi _2}{4}}(((_r+r^1)w^{})^2+c.c.)f_1|w|^2+|w|^4]`$ (21)
Substituting vortex solution $`w=W(r)\mathrm{exp}(i\varphi _0)`$ in Eq. (21), one obtains after integration (since the calculation of the vortex energy is rather straightforward, we refer interested readers to Ref. , p. 11).

$$U=f_1f_2(\overline{\xi }+\xi _2/2\mathrm{cos}2\varphi _0)\mathrm{log}R/r_0+const$$
(22)
where $`r_0O(1)`$ has the meaning of core radius, and $`R`$ is the outer cutoff radius of integration. As one sees from Eq. (22), for physically realizable case $`\xi _2>0`$, the vortices with $`\varphi _0=\pm \pi /2`$ have lower energy, and therefore are more energetically favorable, and are selected in dynamics. In a general case, the vortex solution would have a radius-dependent phase $`\varphi _0`$. Near the center where the density is almost constant, the phase would be close to $`\pm \pi /2`$, and near the island border where the density decreases rapidly, the phase should approach $`0`$ or $`\pi `$. This scenario suggests that the azimuthal velocity should grow with radius near the core, and decrease near edge of the vortex, which is confirmed by our numerical simulations (see below) and agrees with experiments. In addition, one can show that for $`\beta >0`$ the term $`\beta \rho `$ considered as a small perturbation, leads to the drift of the vortex core towards the gradient of density, and therefore it stabilizes the vortex core near the center of the island.

Numerical simulations of the Cahn-Hilliard equation (8) were performed using an FFT split-step method, and the Ginzburg-Landau equation (13) was solved using explicit method. The domain of integration was $`100\times 100`$ dimensionless units with periodic boundary conditions, number of mesh points/FFT harmonics was $`256\times 256`$. As initial conditions we used $`\rho \delta `$ with with small amplitude noise and random initial conditions for $`\psi `$. Selected results are presented in Fig. 2-5.

At the initial stage of the evolution many vortices and small dense clusters (islands) are created throughout the domain of integration (Fig. 2a). Islands are seen as darker areas on the figure because an increase in density $`\rho `$ results in the decrease of the amplitude of $`\psi `$. Some islands trap vortices and are practically immobile, others don’t contain vortices and drift in the direction defined by the mean value of the tilt $`𝐧`$. With time, small islands disappear and bigger islands grow (Fig. 2b,c). It is interesting to note that due to tilt-driven drift coarsening occurs much faster than in ordinary Cahn-Hilliard dynamics. Finally, one big island with the vortex in the center is formed (Fig. 2d). Corresponding density $`\rho `$ and the phase field $`\mathrm{arg}\psi `$ are shown in Fig. 3. Surprisingly, even far from the island there some “dormant vortices” in the low-density phase (gas). These vortices are seen as the end-points of the phase singularity lines (lines between dark and white) in Fig. 3b. These vortices do not annihilate because in the low-density phase the diffusion terms in Eq. (13) are absent. They are not seen in Fig. 3c because their core size is small in the gas phase. However, these vortices can be advected into the high-density phase from the edges.

The velocity field is calculated using Eq. (3). Fig. 4 shows the velocity field corresponding to the calculated tilt structure of Fig.3. The rods perform circular motion around the vortex center. The azimuthal velocity $`v_\theta `$ vs $`r`$ is shown in Fig. 5. The velocity is maximal somewhere between the core and the island edge. It qualitatively resembles the experimental one . Outside the island the tilt becomes large and the velocity becomes small.

In conclusion, we developed a phenomenological model of the formation of the vortical ordered state in the system of vertically vibrated rods. Our continuum model is based on a Cahn-Hilliard equation for the rods aerial density coupled to the Ginzburg-Landau equation for the rod tilt. The model reproduces the qualitative features of the vortex formation process observed in the recent experiment . We thank Toni Neicu and Arshad Kudrolli for stimulating discussions. This work was supported by the U.S.Department of Energy under grants W-31-109-ENG-38 and DE-FG03-95ER14516. Simulations were performed at the National Energy Research Scientific Computing Center.