An optimization parameter for seriation of noisy data
Jeannette Janssen Mahya Ghandehari




























[
    Sandro Mattarei
    December 30, 2023
=====================





 A square symmetric matrix is a Robinson similarity matrix if entries in its rows and columns are non-decreasing when moving towards the diagonal. A Robinson similarity matrix can be viewed as the affinity matrix between objects arranged in linear order, where objects closer together have higher affinity. We define a new parameter, , which measures how badly a given matrix fails to be Robinson similarity. Namely, a matrix is Robinson similarity precisely when its  attains zero, and a matrix with small  is close (in the normalized ℓ^1-norm) to a Robinson similarity matrix. Moreover, both   and the Robinson similarity approximation can be computed in polynomial time. Thus, our parameter recognizes Robinson similarity matrices which are perturbed by noise, and can therefore be a useful tool in the problem of seriation of noisy data. 




 Robinson similarity matrices, Robinsonian matrices, unit interval graphs, seriation, linear embeddings of graphs




  68R01, 05C85, 05C50




§ INTRODUCTION

Many real-life networks, such as online social networks, biological networks and neural networks, are manifestations of an underlying (hidden) spatial reality. For example, members of a social network can be identified with nodes placed in a metric space, whose coordinates represent the interests, backgrounds, and other significant features of the users. The formation of the network is then modeled as
a stochastic process, where the probability of a link occurring between two nodes decreases as their
metric distance increases. A fundamental and challenging problem in the analysis of a social network (or any
other spatial network) is to uncover its “hidden spatial layout”, i.e. to identify the metric space
representation of the network. Analysis and visualization of data becomes considerably more tractable when the dataset is presented according to its spatial reality.

The classical seriation problem, introduced by Robinson in <cit.>, can be viewed as the special case of the spatial layout problem, restricted to one dimension.
The objective of the seriation problem is to order a set of items so that similar items are placed closer to each other. The seriation question translates in a natural way into a question regarding symmetric matrices.
A symmetric matrix is a Robinson similarity  matrix, or Robinson matrix for short, if its entries are non-decreasing when moving towards the main diagonal in each row or column. A symmetric matrix A is said to be  Robinsonian, if it becomes a Robinson matrix after simultaneous application of a permutation π to its rows and columns. In that case, the permutation π is called a Robinson ordering  of A. If the entries of the symmetric matrix A=[A_i,j] represent  similarity of items i and j,  then the  Robinson ordering represents a linear arrangement of the items so that similar items are placed closer together.

The problem of recognizing  Robinsonian  matrices, and finding their Robinson orderings, can be solved in polynomial time.  See <cit.> for the first polynomial time algorithm for this problem, and <cit.> for more recent efficient algorithms. Most of these algorithms are based on a similar principle; namely the connection between Robinsonian similarity matrices and unit interval graphs (<cit.>) or interval (hyper) graphs (<cit.>). A spectral algorithm based on reordering the matrix according to the components of the second eigenvector of the Laplacian, or the Fiedler vector, was given in <cit.>, and was then applied to the ranking problem in <cit.>.

The seriation problem has diverse and significant applications, from its origin in archeological studies to recent applications to ecology and sociology. For a historical overview of the problem and its diverse applications, see <cit.>. 
In most of these applications, it is natural to expect the data to be noisy. In that case, the optimal reordering of a data-derived matrix may not be itself a Robinson matrix, but it will be close to one. The question then becomes, to which extent a given matrix resembles a Robinson matrix. 
All of  the algorithms mentioned in the previous paragraph only apply  to noise-free Robinsonian similarity matrices. 

In the presence of noise, the goal of the seriation problem is to find an “almost Robinson” ordering of a given matrix, i.e. an ordering for which the reordered matrix is closest to being Robinson. This question turned out to be much more challenging than the error-free analogue. In fact,  it is shown in <cit.> that the problem of finding a reordering and a Robinson matrix which is the best ℓ^∞-approximation  is NP-hard. In <cit.>  a factor 16  approximation algorithm is given for the case of ℓ^∞. NP-hardness for a number of related problems is established in <cit.>, where approximation by ℓ^p distance is considered. Specifically, it is shown that for an integer p, the problem of finding proper strong Robinson relations within specified ℓ^p distances of a given matrix is NP-complete. A proper strong Robinson relation corresponds to an appropriate relabelling of the original matrix together with a Robinson matrix with certain additional, stronger properties. 

In <cit.> (together with the references therein) a statistical approach to the problem of seriation with noise is developed, where the error of  Robinson approximation is measured by the Frobenius norm.

If the appropriate labelling is given, then the problem becomes more tractable. Now the problem is that of finding a Robinson matrix closest to a given matrix. For the ℓ^∞-norm, it is known that this problem can be solved in polynomial time since an explicit closed form for the
optimal solution can be easily given (cf. <cit.>). The problem of finding the best ℓ^1-approximation can be formulated as a linear program: Minimize the linear function A-R_1 subject to the constraint that R is Robinson similarity. The constraint can be expressed with O (n^3) inequalities, and thus the problem can be solved in polynomial time.

In this article, we develop new methods and algorithms which can be used for seriation of noisy data. Our focus here is on formalizing the notion of a matrix being “almost Robinson.” To do so,  we introduce a parameter, which we call , that measures how much the local structure of a matrix resembles being a Robinson similarity. Namely,   sums the magnitude  of local violations to the Robinson similarity property, and achieves the value 0 precisely when the matrix is Robinson.  This parameter is a natural tool for the seriation problem, and has been used in practice as a heuristic to measure the amount
of deviation from a Robinson form (see <cit.>).
Moreover,   is simple to formulate, and can be computed in linear time. The main goal of this paper is to show that if the number and magnitude of local violations is small, then the matrix is indeed close, in the sense of  ℓ^1-norm, to  a Robinson matrix.  Precisely, we prove in Theorem <ref> that for every given matrix A,  there exists a Robinson matrix R so that A-R_1≤26(A)^1/3.
In addition, we give a polynomial time algorithm to compute the Robinson approximation R which fulfills the above inequality. 

A proposed application of this work is a novel scheme for treating the seriation problem of noisy data, in which the selection of the best permutation is guided by parameter . The traditional formulation is: Given a matrix A, find a permutation  π and a Robinson matrix R so that A^π-R _p is minimized. (Here A^π refers to the matrix obtained by permuting rows and columns of A according to π.) This approach requires the simultaneous optimization of both the matrix R and the permutation π.
Using the results in this paper, we can instead reduce the noisy seriation problem to the following problem:

Given a symmetric matrix A, find a permutation of its rows and columns so that (A) is minimized. 

Once such a permutation is found, our algorithm can be used to compute the appropriate Robinson approximation. While this approach is not an approximation algorithm per se, we do bound the performance of this approach in terms of the optimal outcome. Namely, as will be shown in Lemma <ref>, the best possible Robinson approximation has normalized ℓ_1 distance at least 1/4 (A) from A.  Therefore, our results implicitly bound the Robinson approximation achieved by our algorithm in terms of the optimal solution. 


The results of this article are fundamentally different from any previous results on ℓ^∞-fitting Robinsonian structures (see for example  <cit.>). Indeed, when matrices grow large in size, the ℓ^1-norm provides us with a more suitable notion of “closeness”.  This fact becomes apparent when we analyze a growing sequence of graphs which are convergent in the sense of Lovász-Szegedy <cit.>. Indeed, this article (and the  choice of notation for the parameter ) was motivated by our previous work <cit.>, where we introduce a parameter Γ which characterizes Robinson graphons. Graphons are symmetric functions on [0,1]^2 with values in [0,1], which can be thought of as the “blueprint” of a random graph whose vertices are randomly sampled from the interval [0,1]. A matrix A=[A_i,j] can be interpreted as a graphon in a natural way, by splitting the unit square  into  subsquares of size 1/n×1/n, and setting the graphon equal to A_i,j everywhere in the (i,j)-th subsquare. A Robinson graphon is a graphon  which is non-decreasing along every horizontal or vertical line towards the main diagonal. 

Our main result in <cit.> is that Γ becomes a continuous parameter, when the space of graphons is equipped with the box-norm. Therefore Γ provides us with a parameter to measure Robinson resemblance, which can be efficiently approximated. The parameters Γ_1 and Γ are closely related, even though box-norm continuity does  not hold for Γ_1 anymore. In future work <cit.>, we employ   these parameters simultaneously in order to develop (continuous) methods for seriation of noisy data. 




Finally, we mention applications to graphs and networks. 
Binary Robinson matrices correspond to unit interval graphs. More precisely, a graph is a unit interval graph if and only if the adjacency matrix is Robinsonian, that is, there exists a labelling of the vertices so that the resulting adjacency matrix is a Robinson matrix. The parameter  can be directly applied to a (labelled) graph, and if  is sufficiently small, our algorithm constructs  a unit interval graph that is close, in edit distance, to the original  graph. 


When applied to real-life networks, the parameter  can be used to measure how closely the matrix conforms to a linear model. In previous work, the authors investigated this question in the context of graph limits <cit.>. In a  linear graph model the vertices of the graph are placed on a line, and the links are formed stochastically so that vertices that are closer together are more likely to connect. The linear layout can refer to a time line, in case of graphs derived from archeology, or the food chain, in case of food webs. But a linear layout may also point to the presence of a strong hidden variable that influences link formation, such as hierarchy, in a professional social network, or age in a friendship graph. If a graph conforms to a linear model,  then we expect the adjacency matrix of the graph to be “almost-Robinson.” Parameter  may therefore serve as a measure of the “linearity" of a given network.   

The rest of this article is organized as follows. In Section <ref>, we introduce the necessary notations and definitions, and we state our main result, namely Theorem <ref>. In Section <ref>, we present Algorithm <ref> which finds a Robinson approximation for the special case where A is a binary matrix (i.e. a graph adjacency matrix). The values of every cell in the Robinson approximation is decided based on the entries of A in the upper right and the lower left regions defined by that cell (see Figure <ref>).  Algorithm <ref> is very simple to state, however one needs to use careful approximations and counting tricks to prove that the algorithm generates an output which indeed is a good ℓ^1-approximation for the input matrix.  Section <ref> provides us with an adaptation of Algorithm <ref> to general matrices. This is indeed a natural generalization, as every matrix with entries in [0,1] decomposes into a convex combination of binary matrices. Fortunately, the parameter  distributes over such decompositions, even though it is not a linear parameter in general. This allows us to apply Algorithm <ref> to the components of the decomposition in stages.  

In Section <ref>, we develop an algorithm that represents a preprocessing step for Algorithm <ref>. This preprocessing step is designed to transform the input matrix A, so that  Algorithm  <ref> generates a better-approximating output matrix R.
The trade-off here is that the preprocessing step increases the complexity of the algorithm, which still remains polynomial. 
We finish the paper by some concluding remarks and future directions. 





§ DEFINITIONS AND MAIN RESULTS

For a given positive integer n, let 𝒜_n denote the set of all symmetric n× n matrices with entries in [0,1].
Note that the restriction on the range of the entries is not a limitation, since it can always be achieved by shifting and scaling of the matrix. 
A matrix A∈𝒜_n is called binary if it has entries only from {0,1}. We refer to the position in the i-th row and j-th column of A as the (i,j)'th cell. For a matrix A of size n, we define its (normalized) ℓ^1-norm to be 
A_1=1/n^2∑_i,j=1^n |A_i,j|. 


An n× n symmetric  matrix A is a Robinson matrix if, for all 1≤ i<j<k≤ n, 

    A_i,j≥ A_i,kA_j,k≥ A_i,k.



In this section, we define a parameter, denoted by , which measures how badly a matrix fails to be Robinson. 
The choice of notation for  is due to the fact that it simply adds the magnitude of violations to (<ref>). 
Precisely, given a symmetric matrix A of size n,

    (A) = 1/n^3∑_1≤ i<k<j≤ n [A_i,j-A_i,k]_+ +[A_i,j-A_k,j]_+,

where [x]_+=x if x≥ 0, and 0 otherwise.

It is clear that (A)=0 if and only if A is a Robinson matrix. Note also that, for a binary matrix A which is not Robinson, (A)≥1/n^3. Moreover,  the computation 
of (A) for an n× n matrix A involves a simple summation which can be executed in O(n^3) steps. Finally note that, due to the normalization factor 1/n^3, (A)∈ [0,1) whenever A∈𝒜_n.

We are now ready to state our main result (Theorem <ref>), whose proof takes a large part of the paper and finishes only at the end of Section <ref>.



For every A∈𝒜_n, there exists a Robinson matrix R∈𝒜_n so that

    A-R_1≤26 (A)^1/3.

Moreover, R can be computed in polynomial time. In addition, if A is binary, then
there exists a binary matrix R satisfying the conditions of the theorem.



The distance between A and the Robinson  approximation obtained in Theorem <ref> is bounded in terms of (A), and thus the algorithm does not necessarily give the best possible approximation. However, as stated in the following simple lemma, there is a close relationship between (A) and the distance between A and the best possible Robinson approximation. 


For every pair of matrices A and R in 𝒜_n, if R is a Robinson matrix then

    A-R_1≥1/4 (A).

Consequently, we have


    min{ A-R_1: R∈𝒜_n }≥1/4 (A).




Since the function [ ·  ]_+ is sub-additive, we observe that 


    (A)=(A-R+R)≤(A-R)+(R)=(A-R),


as R is a Robinson matrix and thus (R)=0. Moreover, for every symmetric matrix B of size n, we have 

    (B)    =   1/n^3∑_1≤ i<k<j≤ n [B_i,j-B_i,k]_+ +[B_i,j-B_k,j]_+
       ≤   1/n^3∑_1≤ i<k<j≤ n |B_i,j|+|B_i,k|+|B_i,j|+|B_k,j|
       ≤   4/n^3∑_k=1^n∑_1≤ i,j≤ n |B_i,j|≤ 4B_1.

Letting B=A-R finishes the proof. 



Applying Theorem <ref> to binary matrices, we obtain an interesting corollary for graphs. For a graph G, let the augmented adjacency matrix B_G denote the matrix which is obtained from the adjacency matrix of G by replacing its diagonal entries by 1. Then, B_G is a Robinson matrix precisely when  the 1-entries of each row and each column are all consecutive. This is called the symmetric consecutive ones property (i.e. C1P for rows and columns simultaneously), and it is known to  characterize unit interval graphs, or equivalently as shown in <cit.>, proper interval graphs (see <cit.>). Precisely, a graph is a unit interval graph if and only if there exists a linear order on its vertices with respect to which B_G has the consecutive ones property, i.e. is a Robinson matrix. 

The parameter  of the augmented adjacency matrix counts the number of triples (i,j,k), 1≤ i<j<k≤ n, for which  vertices i and k are adjacent, but vertex j is not adjacent to either i or k. On the other hand, the (unnormalized) ℓ^1-distance B_G-B_G_1 between the augmented adjacency matrices B_G and B_G of two labeled graphs G and G of the same order corresponds to the edit distance between the graphs themselves. The  edit distance between two graphs G and G, denoted by ed(G,G), is the minimum number of edge deletions and edge additions that need to be performed on G to transfer it to G. Applying these concepts, Theorem <ref> directly leads to the following corollary.



For every graph G  on vertex set V={ 1,2,… ,n}, there exists a unit interval graph G on vertex set V so that 

    ed(G, G)/n^2≤26 (G)^1/3.






§ ROBINSON SIMILARITY APPROXIMATION FOR BINARY MATRICES


In this section, we present an algorithm that finds a Robinson approximation for the special case where A is a binary matrix, and can thus be
 interpreted as the adjacency matrix of a graph. The algorithm can be intuitively understood as follows. We divide all cells of the matrix into black and white cells, and convert all zeros in the black cells to ones, and all ones in the white cells to zeros. The black region is convex around the diagonal, in the sense that, if a cell is black, then so are all other cells closer to the diagonal in the same row or column. In other words, the binary matrix whose support is precisely the black region is a Robinson matrix, which is indeed the Robinson approximation that the algorithm returns (See Figure <ref>).




The decision on whether to assign a cell to the black or white region depends on the entries in the upper right (UR) and lower left (LL) regions defined by the cell. 

Precisely, for any cell (a,b), 1≤ a<b≤ n,  we define

    (a,b)   ={ (i,j):i<a<b<j},
    (a,b)   ={ (i,j):a≤ i≤ j≤ b}.

Roughly speaking, a cell will be black if it has enough ones in its upper right region, and it is white when it has 
enough zeros in its lower left region. So, we need the following notations (also shown in Figure <ref>):

    1_(a,b)   =   |(a,b)∩{ (i,j):A_ij=1}|,
    
    0_(a,b)   =   |(a,b)∩{ (i,j):A_ij=0}|.

In addition, define 1_(a,b)=0_(a,b)=0 when a∈{0,n+1} or b∈{0,n+1} or a>b. 








Algorithm <ref> has complexity θ (n^2), and is therefore linear in the size of the input. A cell is considered black if R_i,j is set to one, and white if R_i,j is set to zero. If (i,j) is black and 
i≤ k<j, then (i,j)⊂ (i,k), so 1_ (i,k)≥ 1_ (i,j)≥ t, and thus (i,k) is also black. Similarly, (k,j) is also black. So, the region of black cells is convex around the diagonal, and R is indeed a Robinson matrix. We will now show in Theorem <ref> that the distance between R and A is bounded by a function of t and n, if the parameter t satisfies Condition (<ref>). 
We will then show in Corollary <ref> that an appropriate t can be chosen as a function 
of (A), so that A-R_1 can be bounded in terms of (A).




Let A∈ be a binary matrix, with the property that 

    for all 1≤ i≤ j≤ n,  1_(i,j) < t or 0_(i,j)< t.

If R is the matrix produced as output of Algorithm <ref> on input A  with threshold t, then A-R_1≤16√(t)+4/n.



As explained earlier, black cells are exactly the cells (i,j) for which R_i,j attains  1, i.e. the cells for which 1_(i,j)≥ t.  Let  denote the collection of all black cells above the diagonal, that is 

    ={(i,j): 1≤ i≤ j≤ n  1_(i,j)≥ t}.


Note that  is convex around the diagonal, in the sense that if a cell (i,j) belongs to  then (i,j)⊆. However, the black region can be disconnected. Precisely, there can be diagonal cells not in , in which case  consists of a collection of connected regions which are convex around the diagonal. Note also that, by definition,  cannot contain any cells (i,j) so that |(i,j)|<t, and specifically,  cannot contain any cells from the first row or the last column.

Let ∂ denote the set of boundary cells of , i.e.

    ∂ ()={(i,j)∈: (i-1,j)∉  (i,j+1)∉}.


The set ∂ precisely contains all black cells above the diagonal that are adjacent to cells outside . Thus, 


    = ⋃_(i,j)∈∂  (i,j).


Similarly, the white region  and its boundary cells are defined as 


    ={(i,j): 1≤ i≤ j≤ n  1_(i,j)< t},

and 

    ∂ ={(i,j)∈:  (i+1,j)∉  (i,j-1)∉}.






Claim 1. Let _0={(i,j)∈: A_i,j=0}. Then |_0|≤ 4n√(t)
 
First, list elements of ∂ as (i_1,j_1),(i_2,j_2),…, (i_m, j_m) in such a way that  i_1≤ i_2≤…≤ i_m, and if i_l=i_l+1 then j_l<j_l+1. This ordering follows the “contour" of the black region, starting with the first black cell in the first row containing any black cells. Since both indices range from at least 1 to at most n, it is clear that the boundary contains at most 2n cells, and m≤ 2n.

To control  the amount of possible overlaps in the covering of  in (<ref>), we now  construct a subsequence C of ∂, so that the lower left regions of cells in C cover most of the black region. 

The subsequence  C={ (i_n_k,j_n_k) } is constructed inductively. In the first step, let ( i_n_1,j_n_1) be the last cell in the first row of .  At step k≥ 1, let n_k be the largest index so that (i_n_k,j_n_k)∈. Define n to be the first index in {1,…, m} which satisfies j_n_k<j_n, (i_n,j_n+1)∉, and either i_n-i_n_k> ⌊√(t)⌋  or j_n-j_n_k> ⌊√(t)⌋. So cell (i_n,j_n) is the first cell in ∂ whose row or column index differs by at least ⌊√(t)⌋ +1 from the last cell added to C, and also is the last black cell in its row. Set n_k+1=n, and add (i_n_k+1,j_n_k+1) to . Since ∂ has at most 2n elements, and ⌊√(t)⌋ +1 ≥√(t), this inductive process ends in at most 2n/√(t) steps.  So ||≤2n/√(t).

We now claim that ∖⋃_(i,j)∈(i,j) can be covered with || squares of dimensions ⌊√(t)⌋×⌊√(t)⌋. Consider two consecutive elements of , say (i_n_k,j_n_k) and (i_n_k+1,j_n_k+1). Let (i,j) be the first cell in ∂∖⋃_(i,j)∈(i,j) after (i_n_k,j_n_k).  By construction, (i,j) is not in the same row or column as (i_n_k,j_n_k), and thus i>i_n_k and j>j_n_k. Let (i',j') be the last cell before (i_n_k+1,j_n_k+1) in ∂∖⋃_(i,j)∈(i,j). (If no such n exists, then k is the last index in C, and we let (i',j')=(i_m,j_m), the last cell of ∂.)


By construction of 
, we know that i'-i+1≤ i'-i_n_k≤⌊√(t)⌋ and j'-j+1≤ j'-j_n_k≤⌊√(t)⌋. Moreover, since (i,j) and (i',j') belong to ∂, and  is convex around the diagonal, every cell (a,b)∈
with  i≤ a≤ i' must satisfy j≤ b≤ j',  so it must belong to the 

square whose top-left corner is (i,j), and its bottom-right corner is (i',j'). We denote this square by S_k, and note that | S_k| ≤ t.

Thus,

    ⊆⋃_k=1^|C|(i_n_k,j_n_k)∪⋃_k=1^|C| S_k.


By Condition (<ref>), every cell (i,j)∈ satisfies 0_(i,j)≤ t. So,


    |_0|≤∑_(i,j)∈ 0_(i,j)+∑_(i,j)∈ |S_i,j|
     ≤ |  | t + |  | t= 4n√(t).





Claim 2. Let _1={(i,j)∈:  A_i,j=1}. Then |_1|≤  4n√(t)+2n.

This proof is similar to the proof of the previous claim. Cells in the boundary ∂ are enumerated similarly, and  a subsequence ={(i_n_k,j_n_k)} of ∂ is defined analogous to : let ( i_n_1,j_n_1) be the last cell in the first column of , and for k≥ 1, define n_k+1 to be the first index in {1,…, m} which satisfies  i_n_k+1>i_n_k, (i_n_k+1+1,j_n_k+1)∉, and  either i_n_k+1-i_n_k> ⌊√(t)⌋ or j_n_k+1-j_n_k> ⌊√(t)⌋. Let S_k be the square of size at most ⌊√(t)⌋×⌊√(t)⌋ covering all white cells between consecutive cells of D. Then we have



    ⊆∂∪⋃_(i,j)∈(i,j)∪⋃_k=1^|D| S_k


Note that  dealing with the white region is slightly different from the black region, in the sense that (i,j) does not include cell (i,j) and the cells in row i and column j (see Figure <ref>). So,  we include the boundary explicitly to cover the white region.




By definition of , every white cell  (i,j)∈ satisfies 1_(i,j)<t. This implies that 


    |_1|≤  |∂ |+ ∑_(i,j)∈ 1_(i,j)+∑_i=1^|| |S_k|≤ 2n+| | t+ || t≤ 4n√(t)+2n.




The above two claims show that R and A differ above the diagonal in at most 8n√(t) +2n cells. Adding the region below the diagonal and normalizing, we conclude that

A-R_1≤2(8√(t)+2)/n.


The following simple counting lemma provides a threshold, in terms of , for which Condition (<ref>) always holds. 


Let A∈ A_n be a binary matrix whose diagonal entries are all 1. 
Then, for every cell (i,j),
we have 

    1_(i,j)  0_(i,j)≤ 2n^4(A).




Without loss of generality, assume that i<j.

    1_(i,j)  0_(i,j)   =   ∑_[   (s,t)∈ (i,j); (s',t')∈ (i,j) ][A_s,t-A_s',t']_+
       ≤   1/2∑_
    1≤ s< i≤ s'≤ t'≤ j< t ≤ n
      [A_s,t-A_s,t']_++[A_s,t'-A_s',t']_+ + [A_s,t-A_s',t]_++[A_s',t-A_s',t']_+
       ≤   1/2∑_
    1≤ s< s'≤ t'< t ≤ n
      [A_s,t-A_s,t']_++ [A_s,t-A_s',t]_+
       +   1/2∑_
    1≤ s<s'≤ t'< t ≤ n
     [A_s,t'-A_s',t']_+ +[A_s',t-A_s',t']_+,


where the second inequality can be justified using the fact that if  [A_s,t-A_s',t']_+=1 then both [A_s,t-A_s,t']_++[A_s,t'-A_s',t']_+=1 and  [A_s,t-A_s',t]_++[A_s',t-A_s',t']_+=1 are satisfied.
Now observe that 

    ∑_
    1≤ s< s'≤ t'< t ≤ n
      [A_s,t-A_s,t']_+≤∑_s'=1^n∑_
    1≤ s< t'< t ≤ n
      [A_s,t-A_s,t']_+≤ n^4(A),


and similarly ∑_
1≤ s< s'≤ t'< t ≤ n
  [A_s,t-A_s',t]_+≤ n^4(A). Moreover, for every s and s', we have [A_s,s'-A_s',s']_+=0, since A_s',s'=1. Thus, 


    ∑_
    1≤ s<s'≤ t'< t ≤ n
     [A_s,t'-A_s',t']_+   =   ∑_
    1≤ s<s'< t'< t ≤ n
     [A_s,t'-A_s',t']_+
       ≤   ∑_t=1^n∑_
    1≤ s<s'< t'≤ n
     [A_s,t'-A_s',t']_+≤ n^4(A),


and similarly ∑_
1≤ s<s'≤ t'< t ≤ n
[A_s',t-A_s',t']_+≤ n^4(A). This finishes the proof. 



For every binary matrix A∈𝒜_n, there exists a binary Robinson matrix R∈𝒜_n such that 

    A-R_1≤5/n+2^9/2(A)^1/4.

Moreover, R can be computed in linear time.


From the binary matrix A∈𝒜_n,  we first construct A by replacing every 0 entry on the diagonal of A with 1. Clearly, A∈𝒜_n and (A)≤(A). Moreover, A-A_1≤1/n. If (A)=0, we take R=A and we are done. So assume that (A)>0. 
Let   R∈𝒜_n be the matrix produced as output of Algorithm <ref> on input A  with threshold t=n^2√(4(A)).
By Lemma <ref> and the fact that (A)>0, A and parameter t=n^2√(4(A)) satisfy Condition (<ref>) of Theorem 
<ref>. Thus,

    R-A_1≤A-A_1+A-R_1≤5/n+16(4(A))^1/4≤5/n+2^9/2(A)^1/4.
 




§ ROBINSON SIMILARITY APPROXIMATIONS OF GENERAL MATRICES

For general matrices, we first decompose the matrix into a convex combination of binary matrices, a standard technique widely used in the literature, and then apply Algorithm <ref> to each summand. Given any matrix A∈_n, let range(A)= { A_i,j:1≤ i≤ j≤ n}. Consider the linear ordering 0=s_0<s_1<… <s_m of  range(A)∪{ 0}, and define matrices A^(k), 1≤ k≤ m as follows:


    A_i,j^(k)=
    {[               1 if A_i,j≥ s_k ,;               0      otherwise. ].

Clearly, each matrix A^(k) is binary, and

    A=∑_k =1^m (s_k-s_k-1)A^(k ).


We refer to the matrices A^(k) as the layers of A. 

From the definition of the layers, it is easy to see that A^(k)≤ A^(l) whenever l<k. Indeed for such l and k, if A^(k)_i,j=1 then 
 A^(l)_i,j=1 as well, since s_l<s_k. 

 
The advantage of writing A as a convex combination of its layers, as opposed to any other decomposition of A, lies in the fact that Γ_1 distributes over this particular decomposition of A, even though Γ_1 is not a linear map in general. 


Suppose A is “layered” as in Equation (<ref>). Then we have

    (A)= ∑_l=1^m (s_l-s_l-1)  (A^(l)).




Fix a triple i,j,k satisfying 1≤ i<k<j≤ n, and let n_1,n_2,n_3∈{1,…,m} be so that 
A_i,j=s_n_1, A_i,k=s_n_2, and A_k,j=s_n_3. From the definition of the layers in (<ref>), we have 


  (i) A_i,j^(l)=1 if l≤ n_1, and A_i,j^(l)=0 if l> n_1.

  (ii) A_i,k^(l)=1 if l≤ n_2, and A_i,k^(l)=0 if l> n_2.

  (iii) A_k,j^(l)=1 if l≤ n_3, and A_k,j^(l)=0 if l> n_3.


Note first that [A_i,j-A_i,k]_+=[s_n_1-s_n_2]_+=0  precisely when n_1≤ n_2. Using (i), (ii) and (iii) it is easy to see that, if n_1≤ n_2, then [A_i,j^(l)-A_i,k^(l)]_+=0 for every 1≤ l≤ m. On the other hand, if [A_i,j-A_i,k]_+>0 and thus n_1> n_2, then, if l≤ n_2 or l>n_1 then  [A_i,j^(l)-A_i,k^(l)]_+=0, and if n_2<l≤ n_1 then [A_i,j^(l)-A_i,k^(l)]_+=1. 
Hence, we can verify the following claim:


    [A_i,j-A_i,k]_+=∑_l=1^m (s_l-s_l-1)[A_i,j^(l)-A_i,k^(l)]_+.


Indeed, if n_1≤ n_2 then both sides of the above equation are equal to 0. For the case where n_1>n_2, we have


    ∑_l=1^m(s_l-s_l-1)[A_i,j^(l)-A_i,k^(l)]_+
    
    =∑_l=n_2+1^n_1s_l-s_l-1=s_n_1-s_n_2=[A_i,j-A_i,k]_+.


Repeating the above argument, we obtain a similar claim for [A_i,j^(l)-A_k,j^(l)]_+, and consequently we get 


    [A_i,j-A_i,k]_++[A_i,j-A_k,j]_+=∑_l=1^m (s_l-s_l-1)([A_i,j^(l)-A_i,k^(l)]_++[A_i,j^(l)-A_k,j^(l)]_+).


So we have  


    ∑_l=1^m (s_l-s_l-1)  (A^(l))   =   ∑_l=1^m (s_l-s_l-1)(1/n^3∑_1≤ i<k<j≤ n [A^(l)_i,j-A^(l)_i,k]_+ +[A^(l)_i,j-A^(l)_k,j]_+)
       =   1/n^3∑_1≤ i<k<j≤ n∑_l=1^m (s_l-s_l-1)([A^(l)_i,j-A^(l)_i,k]_+ +[A^(l)_i,j-A^(l)_k,j]_+)
       =   1/n^3∑_1≤ i<k<j≤ n [A_i,j-A_i,k]_+ +[A_i,j-A_k,j]_+=(A),

which finishes the proof. 






Algorithm <ref> is essentially a simultaneous execution of Algorithm <ref> for each binary matrix A^(k). 
The quantities 1_^k(i,j) thus refer to the number of ones in the region (i,j) in A^(k). The algorithm has complexity O(n^2m), where m is the size of range(A), and thus m≤ n^2. Since every matrix R^(k) is Robinson similarity, so is their linear combination R. In the following theorem, we will show that there exist thresholds t_1,… ,t_m so that the difference between A and R is bounded by C  (A)^1/4 for some constant C. To avoid anomalous behavior, we assume that our input matrix is not Robinson similarity, i.e. (A)>0.



Let A∈_n, and  (A)>0.
If R is the matrix produced as output of Algorithm <ref> on input A  with thresholds t_k=√(4(A^(k)))n^2 for k=1,… ,m, then 

    A-R_1≤2^9/2(A)^1/4(1+O(n^-1/4)).




First suppose A is a binary matrix. Then by Corollary <ref> applied to A with threshold √(4(A))n^2,  we get A-R_1≤5/n+2^9/2(A)^1/4. 



Next, assume that A∈_n is a general matrix, and let A=∑_k =1^m (s_k-s_k-1)A^(k ) be the decomposition of A into layers of binary matrices as described in Equation (<ref>). Then,  we have (A)=∑_k=1^m (s_k-s_k-1)(A^(k)). To avoid clutter of notation, let ϵ=(A) and  ϵ_k:=(A^(k)). For every 1≤ k≤ m, we apply Corollary <ref> to A^(k) with threshold  
t_k=(4ϵ_k)^1/2n^2 to obtain a Robinson similarity matrix R^(k) such that


    A^(k)-R^(k)_1≤5/n+2^9/2ϵ_k^1/4.
 

So, Algorithm <ref> computes R=∑_k=1^m (s_k-s_k-1) R^(k), which is Robinson similarity as well. Moreover, 


    A-R_1   ≤   ∑_k=1^m(s_k-s_k-1)A^(k)-R^(k)_1
       ≤   2^9/2∑_k=1^m(s_k-s_k-1)ϵ_k^1/4+∑_k=1^m(s_k-s_k-1)5/n
       ≤   2^9/2 ((A))^1/4+5/n,


where in the last inequality we used the fact that the function f(x)=x^1/4 is concave. 

By definition, if A∈_n is not Robinson similarity, then (A)≥1/n^3, and thus n(G)^1/4≥ n^1/4. This completes the proof. 





§ IMPROVEMENT THROUGH PREPROCESSING


Finally, we give an algorithm that represents a preprocessing step for Algorithm <ref>. By Theorem <ref>, the Robinson similarity approximation produced by Algorithm  <ref> is at bounded distance from the input matrix A, provided  that Condition (<ref>) holds. By Lemma <ref>, the condition holds for every matrix, if we choose t>√(2)n^2( (A))^1/2 (see proof of Corollary <ref>.) The preprocessing step is designed to transform the input matrix A, such that the condition holds for a smaller value of t. The modified matrix is then used as input to Algorithm  <ref> with the new value of t, which leads to an output matrix R that is closer to the input matrix. 



The trade-off here is that the preprocessing step increases the complexity of the algorithm. This increase is tolerable, as the complexity still remains polynomial. We give an implementation which is quadratic for binary matrices, but we believe that a more sophisticated implementation would lead to an improvement in the complexity. In the following, the algorithm is given in detail for binary matrices only. It can easily be adapted to general matrices, in much the same way that Algorithm <ref> is adapted from Algorithm <ref>. 

First, we introduce some terminology. Suppose that  matrix A is given, and a threshold value t>0 is fixed.  Let Δ denote the collection of all cells in A which are on or above the diagonal. We call a cell (i,j)∈Δ inverted if 1_(i,j)≥ t and 0_(i,j)≥ t. Thus, there exist inverted cells if and only if Condition (<ref>) is violated. To toggle a cell (i,j) is to set the value of all cells in (i,j) equal to zero, and the values of all cells in (i,j) equal to 1. It is easy, yet important, to observe  that toggling a cell can only decrease (A). 



Let A∈ and t>0 be given. Suppose that (i,j)∈Δ is an inverted cell, and let A denote the matrix obtained from A by toggling (i,j). Then for any fixed triple 1≤ s<k<t≤ n, we have


    [A_s,t-A_s,k]_+≤ [A_s,t-A_s,k]_+   [A_s,t-A_k,t]_+≤ [A_s,t-A_k,t]_+.


Consequently, we get (A)≤(A). 




We only prove the first inequality; the second one can be proved similarly. Fix a triple 1≤ s<k<t≤ n. 

First note that, the first inequality in (<ref>) holds trivially whenever (s,t), (s,k)∈Δ∖((i,j)∪(i,j)), as on these cells the matrices A and A are identical. 
For the cases where (s,t)∈(i,j) or (s,k)∈(i,j), we have A(s,t)=0 or A(s,k)=1. Thus, in both cases, we get  [A_s,t-A_s,k]_+=0, and in particular  [A_s,t-A_s,k]_+≤ [A_s,t-A_s,k]_+. Moreover, it is clear from the definition of the region (i,j) that if (s,t)∈(i,j) then (s,k)∈(i,j). Similarly, if (s,k)∈(i,j) then (s,t)∈(i,j) as well. Putting all these together, we conclude that the desired inequality holds in all cases. 

 Therefore, 


    (A)    =   1/n^3∑_1≤ s<k<t≤ n [A_s,t-A_s,k]_+ +[A_s,t-A_k,t]_+
       ≤   1/n^3∑_1≤ s<k<t≤ n [A_s,t-A_s,k]_+ +[A_s,t-A_k,t]_+
       =   (A),

which finishes the proof. 




Algorithm <ref> can then be described as follows: each inverted cell is toggled, and after each toggling step, the values of 1_ and 0_ are recalculated for each cell. 



Note that Algorithm <ref> involves computing  1_(a,b) for all 1≤ a≤ i and j≤ b≤ n, and computing 0_(a,b) for all i≤ a ≤ j and i≤ b≤ j. These values must be recomputed in each iteration, since A is being changed. Any cell (i,j) is tested exactly once, to see whether it is inverted, and if it is, it is toggled.  This naive implementation of the algorithm has complexity O(n^4) and is thus quadratic in the size of the input. When adapted to general matrices, the complexity becomes O(mn^4), where m is the number of different values taken by entries of A. Clearly m≤ n^2. Also, if all entries of A are rounded to the nearest multiple of ϵ^A1/3 (recall that ϵ =Γ_1 (A)), then m=ϵ^-1/3, while the error between A and R  is still of the same order.

One may be concerned  that toggling some cells may create new inverted cells, in which case, just considering each cell once would not be sufficient, and the complexity could increase. The following lemma shows this cannot be the case. 


Suppose Algorithm <ref> is applied to a binary matrix A, with threshold t. Then  the output of the algorithm, which we call the modified matrix, satisfies the condition that, for each cell (i,j), 1_(i,j)<t or 0_(i,j)<t. 



It suffices to show that, if an inverted cell (i,j) is toggled, no cell that was not inverted before the toggle can become inverted afterwards. Suppose to the contrary that there exists a cell (k,ℓ) which becomes inverted after the toggle. 
That is, either 1_(k,ℓ) is increased by the toggle, or 0_ (k,ℓ) is increased. Assume without loss of generality  that 1_(k,ℓ) increases after the toggle. 
The toggle sets all entries in (i,j) to zero, and all entries to (i,j) to one. So for 1_(k,ℓ) to increase, (k,ℓ) and (i,j) must intersect. Thus, (k,ℓ)∈ (i,j). Therefore, after the toggle, all entries in (k,ℓ) have been set to one, so 0_(k,ℓ)=0<t. Therefore, (k,ℓ) is not an inverted cell after the toggle, which contradicts our assumption.


The above lemma shows that the modified matrix A returned by Algorithm <ref> satisfies Condition <ref> of Theorem <ref> for the chosen value of t. Thus, this matrix can be used as input for Algorithm <ref> to obtain a Robinson similarity approximation R, and Theorem <ref> can be applied to bound A-R_1. However, to obtain a good approximation we need to make sure we can  also  bound  A-A_1, the distance between the input and output matrices of Algorithm <ref>. The following lemma gives such a bound. 



Let A∈_n be a binary matrix. Let A denote the output of Algorithm <ref> with threshold t. Then A-A_1≤4n^2/t((A)-(A)).



The output matrix  A is obtained from A by consecutive toggling of inverted cells, occurring in lines 5 and 6 of Algorithm <ref>. The condition in line 4 checks whether a cell is inverted. Let A=A^0, A^1,A^2,…, A^m=A denote the matrices in the intermediate steps. We will bound the distance between two consecutive matrices.

Fix s, 0≤ s<m, and assume that A^s is modified because cell (i,j) is found to be inverted, and is then toggled.
So 1_(i,j)≥ t and 0_(i,j)≥ t, where 1_ and 0_ are computed from A^s. So A^s+1 is formed from A^s by adjusting every cell in (i,j) and its counterpart below the diagonal to be 0, and every cell in (i,j) and its counterpart below the diagonal to be 1. 

We can intuitively observe that (A^s) drops by at least 0_LL(i,j)1_UR(i,j)/n^4 when (i,j) is toggled. Note that, normalizing as in the definition of  and applying the same reasoning as in Lemma <ref>, one may be tricked into thinking that (A^s) drops at least by 2/n^30_LL(i,j)1_UR(i,j) when (i,j) is toggled, since a contribution of 2/n^3 is removed for each pair of “bad” cells from (i,j) and (i,j). However, this reasoning does not take overcounting into consideration. We dedicate the following claim to a rigorous proof of this intuitive observation. 



For A, s, and (i,j) as above, we have

    (A^s)-(A^s+1)≥0_(i,j)1_(i,j)/n^4.



First, consider a cell (k,ℓ) from (i,j) containing 1,  and a cell (k',ℓ')  from (i,j), containing 0. 
Similar to the proof of Lemma <ref>, 

we have that 

    2=2[A^s_k,ℓ-A^s_k',ℓ']_+    =    [A^s_k,ℓ-A^s_k',ℓ]_+ +[A^s_k',ℓ-A^s_k',ℓ']_+
          +[A^s_k,ℓ-A^s_k,ℓ']_+ +[A^s_k,ℓ'-A^s_k',ℓ']_+.


By Lemma <ref>, for every triple  1≤ k<l'<l≤ n we have 

    [A^s+1_k,l-A^s+1_k,l']_+≤ [A^s_k,l-A^s_k,l']_+   
    [A^s+1_k,l-A^s+1_l',l]_+≤ [A^s_k,l-A^s_l',l]_+.


This, together with the fact that A^s+1_k,l=0 whenever (k,l)∈(i,j), implies that

    ∑_[ 1≤ k<l'<l≤ n;  (k,l)∈(i,j);     i≤ l'≤ j ] [A^s_k,l-A^s_k,l']_+   =   ∑_[ 1≤ k<l'<l≤ n;  (k,l)∈(i,j);     i≤ l'≤ j ] [A^s_k,l-A^s_k,l']_+-[A^s+1_k,l-A^s+1_k,l']_+
       ≤   ∑_1≤ k<l'<l≤ n [A^s_k,l-A^s_k,l']_+-[A^s+1_k,l-A^s+1_k,l']_+.


Similarly, we have 

    ∑_[ 1≤ k<k'<l≤ n;  (k,l)∈(i,j);     i≤ k'≤ j ] [A^s_k,l-A^s_k',l]_+
       ≤   ∑_1≤ k<k'<l≤ n [A^s_k,l-A^s_k',l]_+-[A^s+1_k,l-A^s+1_k',l]_+.


Adding up the above two inequalities, we get,

    ∑_[ 1≤ k<l'<l≤ n;  (k,l)∈(i,j);     i≤ l'≤ j ] [A^s_k,l-A^s_k,l']_++
    ∑_[ 1≤ k<k'<l≤ n;  (k,l)∈(i,j);     i≤ k'≤ j ] [A^s_k,l-A^s_k',l]_+
     ≤ 
    n^3(Γ(A^s)-Γ(A^s+1)).


Repeating the above argument for elements of (i,j), the following inequality can be derived in a  similar fashion. 


    ∑_[ 1≤ k<k'<l'≤ n; (k',l')∈(i,j);       1≤ k≤ i ] [A^s_k,l'-A^s_k',l']_++
    ∑_[ 1≤ k'<l'<l≤ n; (k',l')∈(i,j);       j≤ l≤ n ] [A^s_k',l-A^s_k',l']_+
     ≤ 
    n^3(Γ(A^s)-Γ(A^s+1)).

Combining (<ref>), (<ref>), and (<ref>), we conclude the following.


    1_(i,j) 0_(i,j)   =   ∑_[   (k,l)∈(i,j); (k',l')∈(i,j);               ]
    [A^s_k,ℓ-A^s_k',ℓ']_+ 
       =   1/2∑_[   (k,l)∈(i,j); (k',l')∈(i,j);               ]([A^s_k,ℓ-A^s_k',ℓ]_+ +[A^s_k,ℓ-A^s_k,ℓ']_+)
       +   1/2∑_[   (k,l)∈(i,j); (k',l')∈(i,j);               ]([A^s_k',ℓ-A^s_k',ℓ']_+
     +[A^s_k,ℓ'-A^s_k',ℓ']_+)
       ≤   n/2∑_[ 1≤ k<l'<l≤ n;  (k,l)∈(i,j);     i≤ l'≤ j ] [A^s_k,l-A^s_k,l']_++
    n/2∑_[ 1≤ k<k'<l≤ n;  (k,l)∈(i,j);     i≤ k'≤ j ] [A^s_k,l-A^s_k',l]_+
       +   n/2∑_[ 1≤ k<k'<l'≤ n; (k',l')∈(i,j);       1≤ k≤ i ] [A^s_k,l'-A^s_k',l']_++
    n/2∑_[ 1≤ k'<l'<l≤ n; (k',l')∈(i,j);       j≤ l≤ n ] [A^s_k',l-A^s_k',l']_+
       ≤    n^4(Γ_1(A^s)-Γ_1(A^s+1)).

Note that the factor n/2 in the above inequalities appear, because every fixed row or column of A^s can have at most n cells from (i,j) or (i,j).
This completes the proof of the claim. 



Finally, it is easy to observe that A^s-A^s+1_1≤2(0_(i,j)+1_UR(i,j))/n^2. Without loss of generality, assume that 1_≥ 0_. Then


    (A^s)-(A^s+1)≥0_LL(i,j)1_UR(i,j)/n^4=( 0_ (i,j)/4n^2)( 41_UR(i,j)/n^2)≥t/4n^2A^s-A^s+1_1.


 Applying the above result, we get

    (A)-(A)   =   ∑_s=0^m ((A^s)-(A^s+1))≥∑_s=0^mt/4n^2A^s-A^s+1_1
       ≥   t/4n^2∑_s=0^m A^s-A^s+1_1=t/4n^2A-A_1,


which finishes the proof.



We now have all the ingredients to prove Theorem <ref>, by combining the above lemma with Theorem <ref>. For this, parameter t must be tuned so that the bounds from the preprocessing step and from Algorithm <ref> give the best possible result; it appears that the best choice is  t=4^-2/3(A)^2/3n^2. (Note that this value is smaller than the value used in Corollary <ref>.) 
For this choice of t,  the distance between the input matrix A and the Robinson similarity matrix returned by Algorithm <ref> when applied to the updated version of matrix A (after being processed  by Algorithm <ref>) is bounded by 26 (A)^1/3. As the exponent on (A) has been decreased from 1/4 to 1/3, the preprocessing step leads to a substantially better Robinson  approximation.



First assume that A∈_n is a binary matrix with (A)=ϵ>0. Let A be the output of Algorithm <ref>, with threshold t=4^-2/3ϵ^2/3n^2. From Lemma <ref>, 

    A-A_1≤4^5/3ϵ^-2/3((A)-(A))≤ 4^5/3ϵ^1/3.
 
From Lemma <ref>, we have that A satisfies Condition (<ref>) of Theorem <ref> for our choice of t.
Let R be the output of  Algorithm <ref> applied to A with parameter t=4^-2/3ϵ^2/3n^2. Then, by Theorem <ref>, we have

    A-R_1≤16√(t)+4/n=4^5/3ϵ^1/3+4/n.


Combining these inequalities, we get A-R_1≤ 2· 4^5/3ϵ^1/3+4/n≤26ϵ^1/3, where we used the fact that (A)=ϵ≥1/n^3 in the last inequality. 

Next, assume that A∈_n is a general matrix, and let A=∑_k =1^m (s_k-s_k-1)A^(k ) be the decomposition of A into layers of binary matrices as described in Equation (<ref>). Let ϵ_k:=(A^(k)), and recall that ϵ=∑_k=1^m (s_k-s_k-1)ϵ_k. For every 1≤ k≤ m, we apply the process described in the above paragraph to  A^(k ) with parameter t_k=4^-2/3ϵ_k^2/3n^2, and we obtain  Robinson matrices R^(k) such that


    A^(k)-R^(k)_1≤26ϵ_k^1/3.
 

Letting R=∑_k =1^m (s_k-s_k-1)R^(k), we have


    A-R_1   ≤   ∑_k=1^m(s_k-s_k-1)A^(k)-R^(k)_1
    ≤26∑_k=1^m(s_k-s_k-1)ϵ_k^1/3≤26ϵ^1/3,


where in the last inequality we used the fact that the function f(x)=x^1/3 is concave. 

Finally,
we observe that the outputs of Algorithm <ref> and Algorithm <ref>, when applied to a binary matrix A, are again binary matrices. Therefore, the Robinson approximation R of Theorem <ref> is binary, when A is a binary matrix.
 



§ CONCLUSIONS AND FURTHER WORK

We defined a parameter  which measures how much a matrix resembles a Robinson similarity matrix. We gave a polynomial time algorithm which takes as input a symmetric matrix A, and finds a Robinson matrix R so that the normalized ℓ^1-distance between A and R is bounded by 26 (A)^1/3. The motivation of our work is the application to the problem of seriation of noisy data. This problem can now be approached by solving instead the optimization problem: given a matrix A, find a permutation π of the rows and columns of A so that (A^π) is minimized. 

Our construction method is based on a combinatorial algorithm that runs in polynomial time. However, this may not give the best possible such Robinson approximation. As remarked in the introduction, the problem of finding the best possible Robinson approximation, with error measured in ℓ_1 norm, can be formulated as a linear program. An open problem is whether there exists a combinatorial algorithm for this task (as there exist for the ℓ_∞ norm.)

In future work, we propose to study this optimization problem, to attempt to find algorithms which solve or approximate the problem, and to determine their complexity. It is well-known that the second eigenvector of the Laplacian of the matrix, also known as the Fiedler vector, is effective in finding the correct permutation in seriation without error. We propose to study the relationship between the Fiedler vector and the parameter  . 

An immediate next step is to test our algorithm on real data, and see whether, in practice, the algorithm outperforms the theoretical bound. As well, a clever implementation of Algorithm <ref> will likely lead to improved efficiency.   



§ ACKNOWLEDGMENTS


The authors thank the anonymous referees, whose suggestions greatly improved the paper.


siamplain

10

atkins
J. E. Atkins, E. G. Boman, and B. Hendrickson, A spectral algorithm
  for seriation and the consecutive ones problem, SIAM J. Comput., 28 (1999),
  pp. 297–310, <https://doi.org/10.1137/S0097539795285771>,
  <http://dx.doi.org/10.1137/S0097539795285771>.

Barthelemy-Brucker
J.-P. Barthélemy and F. Brucker, NP-hard approximation
  problems in overlapping clustering, J. Classif., 18 (2001), pp. 159–183,
  <https://doi.org/10.1007/s00357-001-0014-1>,
  <https://doi.org/10.1007/s00357-001-0014-1>.

Chen2002
C.-H. Chen, Generalized association plots: information visualization
  via iteratively correlated matrices, Statistica Sinica, 12 (2002),
  pp. 7–29.

chepoi-fichet
V. Chepoi, B. Fichet, and M. Seston, Seriation in the presence of
  errors: NP-hardness of ℓ_∞-fitting robinson structures to
  dissimilarity matrices, J. Classification, 26 (2009), pp. 279–296.

Chepoi2011
V. Chepoi and M. Seston, Seriation in the presence of errors: A
  factor 16 approximation algorithm for ℓ_∞-fitting robinson
  structures to distances, Algorithmica, 59 (2011), pp. 521–568.

linearembeddings
H. Chuangpishit, M. Ghandehari, M. Hurshman, J. Janssen, and
  N. Kalyaniwalla, Linear embeddings of graphs and graph limits, Journal
  of Combinatorial Theory, Series B, 113 (2015), pp. 162 – 184,
  <https://doi.org/https://doi.org/10.1016/j.jctb.2015.02.002>,
  <http://www.sciencedirect.com/science/article/pii/S0095895615000179>.

uniform
H. Chuangpishit, M. Ghandehari, and J. Janssen, Uniform linear
  embeddings of graphons, European J. Combin., 61 (2017), pp. 47–68,
  <https://doi.org/10.1016/j.ejc.2016.09.004>,
  <http://dx.doi.org/10.1016/j.ejc.2016.09.004>.

Corneil2004
D. G. Corneil, A simple 3-sweep LBFS algorithm for the recognition
  of unit interval graphs, Discrete Applied Mathematics, 138 (2004), pp. 371
  – 379.

Corneil1995
D. G. Corneil, H. Kim, S. Natarajan, S. Olariu, and A. P. Sprague, Simple linear time recognition of unit interval graphs, Information
  Processing Letters, 55 (1995), pp. 99 – 104.

Rigollet
N. Flammarion, C. Mao, and P. Rigollet, Optimal rates of statistical
  seriation.
arXiv:1607.02435 [math.ST], 2016.

fogel
F. Fogel, A. d'Aspremont, and M. Vojnovic, Spectral ranking using
  seriation, J. Mach. Learn. Res., 17 (2016), pp. Paper No. 88, 45.

GhandJanss2019
M. Ghandehari and J. Janssen, A stable parameter to quantify
  geometric structure of graphons.
Unpublished Manuscript, 2018.

Hahsler2008
M. Hahsler, K. Hornik, and C. Buchta, Getting things in order: An
  introduction to the R package seriation, J. Statistical Software, 25 (3)
  (2008).

laurent2017
M. Laurent and M. Seminaroti, A lex-bfs-based recognition algorithm
  for robinsonian matrices, Discrete Applied Mathematics, 222 (2017), pp. 151
  – 165, <https://doi.org/https://doi.org/10.1016/j.dam.2017.01.027>,
  <http://www.sciencedirect.com/science/article/pii/S0166218X17300641>.

laurent2016
M. Laurent and M. Seminaroti, Similarity-first search: a new
  algorithm with application to robinsonian matrix recognition, SIAM
  J. Discrete Math., 31 (2017), pp. 1765–1800,
  <https://doi.org/https://doi.org/10.1137/16M1056791>.

Liiv2010
I. Liiv, Seriation and matrix reordering methods: An historical
  overview, Stat. Anal. Data Min., 3 (2010), pp. 70–91.

Looges1994
P. J. Looges and S. Olariu, Optimal greedy algorithms for
  indifference graphs, Computers & Mathematics with Applications, 25
  (1993), pp. 15 – 25.

lovasz-szegedy
L. Lovász and B. Szegedy, Limits of dense graph sequences, J.
  Combin. Theory Ser. B, 96 (2006), pp. 933–957,
  <https://doi.org/10.1016/j.jctb.2006.05.002>.

mirkin-rodin
B. G. Mirkin and S. N. Rodin, Graphs and genes, vol. 11 of
  Biomathematics, Springer-Verlag, Berlin, 1984.
Translated from the Russian by H. Lynn Beus.

prea-fortin
P. Préa and D. Fortin, An optimal algorithm to recognize
  Robinsonian dissimilarities, J. Classification, 31 (2014), pp. 351–385,
  <https://doi.org/10.1007/s00357-014-9150-2>,
  <http://dx.doi.org/10.1007/s00357-014-9150-2>.

Roberts1969
F. Roberts, Indifference graphs, in Proof Techniques in Graph
  Theory, F. Harary, ed., Academic Publishers, 1969, pp. 139–146.

robinson
W. S. Robinson, A method for chronologically ordering archeological
  deposits, American Antiquity, 16 (1951), pp. 293–301.

Seminaroti2016
M. Seminaroti, Combinatorial Algorithms for the Seriation Problem,
  PhD thesis, Un. Tilburg, 2016.

seston
M. Seston, A simple algorithm for recognize Robinsonian
  dissimilarities, in COMPSTAT 2008—Proceedings in Computational
  Statistics, Physica-Verlag/Springer, Heidelberg, 2008, pp. 241–248,
  CD–ROM.




