Restoring a list from differences

Question

Given a list of (absolute valued) pair differences ordered and with duplicates removed, how can we recover/reconstruct the list that generated these differences? We do not know anything about the generating list, not the individual values and not its length.

Normally, the objective is to find the shortest possible generating list, but also to find all solutions.

Here is an example: Given the (absolute valued) difference list of $$[1,6,7,8,12,13,14,15,19,20,21,26,27,33,34]$$ we find the following solutions, i.e. lists that could have generated these pair differences:

$[1,2,9,16,22,28,35],[1,8,14,20,27,34,35],[1,2,8,16,22,28,35],[1,8,14,20,28,34,35]$

These solutions are "normalized" in that they always start with $1$ and end with $1+\max(\text{difference list})$.

The difference list was actually generated by the list $[5,12,18,24,31,32,38,39]$. It is "normalized" (all numbers are subtracted by some constant integer) to $[1,8,14,20,27,28,34,35]$ which happens to be one of the found solutions.

I have thought about this for some days now and written:

• a Picat (constraint programming) model with some experiments. The model also handles cases where there are duplicates in the difference lists.

• a MiniZinc model implementing the base version. Here is a MiniZinc-Python program for experimenting with the MiniZinc model.

• a draft describing the problem and some experiments/findings/invariants. Unfortunately, right now it's not as "formal" as I would want it.

The base model that is described below has been implemented in the Picat model (link above). The description uses constraint programming parlance.

• "origin" list (generator): This "origin" list is not really in the model but it is essential to understand the model. The "origin" is the list that generates the list of pair differences that is the input of the model. The origin list is not known by the model and we don't know either its length or any of its elements (integers).

  The objective of the model is to find out what this origin list is, or rather some variant of it since in real use we don't know what the origin list is.

  Let us call this origin list L with length Len. Then the difference list is calculated as follows:

  • Get all the pair differences of L: L2 = [abs([L[I]-L[J]) : I in 1..Len, J in I+1..Len]

  • Sorting the pair difference list L2 and removing the duplicates give the sorted difference list: Diffs = L2.sort_remove_dups.

  Note: This part is used in the model just to generate difference lists to test. The input to the model is the difference list.

• "normalization" of the origin: Each origin list can be normalized to a list which starts with $1$. Let L2 be the sorted version of L and with length LLen. Then the normalized version is [L2[I] - L2[1] + 1 : I in 1..LLen], i.e. subtract all integers in the list L2 with the first value of the list and add $1$ (the start value of the normalized list).

  Note: Instead of $1$ as the start value, we could use $0$ instead (or any other value), but I tend to prefer that the list start with $1$.

The concept of normalization is important since the solution of the model can be seen as a normalized version of the origin.

The input of the model is

• diffs: The difference list that is generated by the (unknown) origin list as described above: it is an ordered list of distinct integers. It has length len and a maximum value (the last value) max_val.

The model:

The model loops through a value of lengths of the solution list, starting with min_len:

• min_len: The minimum length of a possible solution. This value is calculated as $$\min\_\text{len}=\left\lfloor1+\frac{\sqrt{1+8\cdot\text{len}}}2\right\rfloor.$$ The formula is based on the calculation of the number of difference pairs in a given list of length $n$, which is given by $n(n-1)/2$. Thus, given the length of the difference list (len), the minimum length is min_len.

• For finding one of the shortest solutions, we loop over the the possible lengths (n) of the solution list, from the minimum length min_len to some (arbitrary) upper length.

  Note: I have not found a way of calculating this upper value given the inputs of the model. The upper value of the length of the solution is actually the length of the origin list, but in the model we don't know anything about the origin list.

  Note 2: In most cases, the length of the solution is min_length or some value slightly larger than that.

• Each loop of possible length n (from min_length) does the following:

  • creates the list of decision variables x of length n.

  • the minimum value of x is 1 and the upper value is max_val + 1, i.e. the domain is $\{1,\cdots,\max\_\text{val}+1\}$. The reason behind this upper value is that if the solution starts with $1$, then the largest possible difference must be the largest difference value ${}+ 1$.

  The solution is thus thought of as a "normalized" version of the origin list, or rather a possible variant since there might be many different solutions where only one of then might be the same as the normalized version of the origin list.

  • add constraints of the list x as described below.

  • if we find a solution then it is a shortest solution. Normally the model is finished after this. But there are cases where we continue to find all the solutions or perhaps all the shortest solutions. (The Picat model does quite a few experiments with this.)

The constraints in the model are:

• all_different(x): Ensures that all the decision variables are distinct. Note that this is not really necessary since the next constraint increasing_strict(x) ensures this as well. For constraint programming solvers, adding this constraint seems to speed up the model.

• increasing_strict(x): Ensures that the decision variables are strictly ordered, i.e. distinct and sorted. The ordering part can be seen as a symmetry breaking to reduce a lot of symmetric solutions. The strict part is necessary for the model though.

  (Note: In some experiments, the strict part is skipped so we can study origin lists and difference lists that contain duplicates. This is not discussed further here.)

• x[1] = 1: the first element is $1$. This is to make the solution "normalized".

• x[n] = max_val: the last element is max_val, i.e. 1 + max(diffs). These two constraints make sure that the solution is "normalized".

• All the pair differences in the solution x must be in the difference list diffs. The pseudo code for this is:

foreach(i in 1..n, j in i+1..n)
   abs(X[I]-X[J]) in diffs
end

i.e. All possible pair differences in x must be in the difference list.

• Ensure that we cover all the differences in the difference list diffs. Pseudo code:

foreach(d in diffs)
   % Create the indices i,j where I < J
   i in 1..n,
   j :: 1..n,
   i < j,
   d = abs(x[i]-x[j])
end

Note that i and j are decision variables with domain $\{1,\cdots,n\}$. For all differences d in the difference list, we try to find some indices i and j (with i < j) such that d = abs(x[i]-x[j]).

The two constraints above ensure that the solution list x generates all the differences in difference list diffs as well as only those differences.

This is the end of the base model. It is - as mentioned above - implemented in Picat and it works quite well. However, it would be interesting to know alternative approaches for solving this problem.

I haven't found much about this, but perhaps the terms I searched for were not the proper ones.

One vaguely related paper is a paper by De Biasi (2014)¹ about reconstruction permutations of length $N$ given a permutation of differences of length $N-1$.

However, it only handles difference lists that are permutations, but I am interested in any difference lists, not just permutations. Also, the objective of the paper is to prove that the problem is NP-complete which is not my goal here. (Here is a Picat model that solves these kind of problems.)

One other related thing is perfect rulers where the difference lists are the full list of integers $1,\cdots,N$:

Perfect rulers might be seen as a special case of what I am interested in, but the general case of difference lists is that they may be non-contiguous. And the objective is not the same: I want to recover the list that generated the differences.

Does anyone have more information about this?

---

_Reference

_{[1] De Biasi, M. (2014). Permutation Reconstruction from Differences. The Electronic Journal of Combinatorics. 21(4):P4.3.}

Answer 1

It will be useful to distinguish the problem itself from the algorithms and data structures needed to solve it. To that end I will begin with a brief mathematical analysis, in which I will use standard mathematical notation without further apology, because all definitions are readily available elsewhere. Non-standard terms specially coined for this analysis will be introduced within quotation marks. Next I discuss a data structure and algorithm for solving the problem with reasonable efficiency. Some remarks about the computational costs follow. This post ends with examples and executable code.

---

Analysis

The difference of two subsets $A,B$ of a group $(G, +)$ is the set of all possible pairwise differences,

$$A - B = \{a-b\mid a\in A,\; b\in B\}.$$

The difference set of any subset $A\subset G$ is just

$$\Delta(A) = A - A.$$

$G$ acts on its subsets $\mathcal{P}(G)$ via translation; the translate of any $A$ by an element $g\in G$ is

$$A + g = \{a+g\mid a\in A\}.$$

Because $\Delta(A+g) = \Delta(A)$ for all $g\in G,$ we may view $\Delta$ as a function defined on the orbits $\mathcal{P}(G)/G,$

$$\Delta: \mathcal{P}(G)/G \to \mathcal{P}(G).$$

Given $X\subset G,$ $[X]\in\mathcal{P}(G)/G$ represents its orbit.

The problem is to invert $\Delta;$ that is,

Given nonempty $D \in \mathcal{P}(G),$ find all orbits $[X]\in \mathcal{P}(G)/G$ for which $\Delta([X]) = D.$

In what follows I will focus on the case $G=(\mathbb{Z},+),$ the integers with their usual addition, but little would need to change to apply the solution to other commutative groups (of which the additive real numbers and multiplicative complex numbers are likely to have interesting applications.)

An Abstract Algorithm

Throughout this section $D$ is fixed and all definitions are relative to $D.$

The proposed approach builds solutions incrementally. To this end, pick once and for all a unique representative $X$ for each orbit $[X].$ (This is the "normalization" discussed in the question.)

Let us say that $X^\prime$ is a "partial solution" to the problem when $\Delta(X^\prime) \subset D.$ Solutions are partial solutions.

In such a case, any $g\in G\setminus X^\prime$ increases $X^\prime$ to a larger partial solution, or "augments" $X^\prime,$ when including $g$ in $X^\prime$ still gives a partial solution. I will call the set of all such $g$ the "augmentors" of $X$ and write it $\mathcal{A}_D(X).$

The import of augmentation is that every solution can be generated as a sequence of augmentations of the empty set $\{\}.$ This is obvious, because if $X$ is a solution it must be nonempty, whence there exists $g\in X$ and clearly (a) $X^\prime=X\setminus\{g\}$ is a partial solution and (b) $X$ augments $X^\prime.$ Thus, by successively removing elements from a solution one by one, ultimately we must arrive at the empty set.

The foregoing defines a graph on the orbits $\mathcal{P}(G)/G.$ Every orbit representative $X\subset G$ has a (possibly empty) set of augmentors $\mathcal{A}_D(X).$ Each $g\in \mathcal{A}_D(X)$ determines a directed edge from $X$ to $X\cup \{g\}.$ Because the cardinality of the nodes increases by $1$ along all such edges, there can be no cycles: this graph is a forest. We may restrict our analysis to the connected component--a tree--rooted at the empty set. All partial solutions, and therefore all solutions, are nodes in this tree. Therefore, the problem is seen to be one of finding all true solutions in this partial-solution tree.

The obvious algorithm for finding all solutions is to traverse the tree.

(Of course if you need only one solution you may stop the traversal upon encountering the first solution.)

The traversal can be made more efficient by remembering which candidates have been tried previously during the traversal. This is readily done by sorting the augmentors and trying them in sort order in a depth-first traversal. After exploring the partial solution $X^\prime$ by following the edge determined by an augmentor $g,$ we may eliminate $g$ from the augmentors of all descendants of $X^\prime$ when we continue the search to other children of $X^\prime.$ The augmentors of the descendants that haven't been eliminated in this fashion are "candidates" in the search for a solution.

Thus, the traversal algorithm depends on the current node $X^\prime$ and a set $C$ of candidates for it. If $X^\prime$ is a solution, output it. To search the subtree rooted at $X^\prime$, for each candidate $g$ recursively apply the algorithm to the node $X = X^\prime\cup\{g\}$ and the candidate set $\mathcal{A}(X)\,\cap\, (C\setminus\{g\}).$

Data Structures

In searching for solutions it will be important to support various operations efficiently, including

1. Representing sets, orbits, and partial solutions.

2. Comparing sets for equality.

3. Adjoining and removing elements from sets.

As a matter of creating efficient data structures, notice

• The elements of $\Delta(X)$ appear in pairs $\{a-b,b-a\}$ and one of these pairs is always $\{0,0\}=\{0\}$ provided $X$ is nonempty.

• In an ordered group like $\mathbb Z,$ the values $\pm(\max(X) - \min(X)) \in \Delta(X)$ are the most extreme elements and therefore must agree with $\pm\max(D)$ when $\Delta(X) = D.$

Accordingly, the basis for a data structure will include an efficient representation of a set $\Delta(X)$ where $X\subset \mathbb{Z}$ is nonempty.

All the necessary information for solving this problem can be maintained in a data structure with the following elements (whose names reflect the origin of this problem in physical spectra, which represent differences between energy levels in a quantized system):

• Delta represents the original putative difference set $D \subset \mathbb{Z}$. It need only include the positive elements of $D.$

• Levels represents a partial solution $X^\prime.$ It needn't begin with the empty set, because any solution must include $0$ and $n.$

• Spectrum represents the difference set of Levels, $D(X^\prime).$

• Candidates represents the candidate set $C.$ It will need to be dynamically updated during the tree traversal.

Example

Consider $D = \{0,1,3,4,5,7,8,9\}.$ We may begin the search at $X^\prime = \{0,9\}$ because $n=9$ is the range of the differences $D$ and I am normalizing the orbit representatives to begin at $0$ (rather than $1$). The initial candidate set $C$ consists of the augmentors of $X^\prime.$ These will be the integers $g$ between $0$ and $9$ (exclusive) for which both $g-0\in D$ and $9-g\in D,$ yielding $C = \{1,4,5,8\}.$ Thus the root node of the search tree is

$$(X^\prime, C) = (\{0,9\},\ \{1,4,5,8\}).$$

In the figure I have dropped the set braces and commas for brevity.

This diagram of the (pruned) tree for $\{0,1,3,4,5,7,8,9\}$ is rooted at the top. The candidates $g$ corresponding to the edges $(X^\prime;*)\to(X;*)$ can be determined by comparing $X$ to $X^\prime:$ $g$ is the new element in $X.$ Partial solution nodes are shown in blue and full solutions in orange.

This is not the full tree: it is the tree after it has been pruned during a depth-first, left-to-right search. Thus, for instance, after initially searching the candidate $g=1$ from the root in the leftmost branch, $1$ was eliminated from the candidate lists in the subsequent branches from the root.

The three solutions are, in the order found, $\{0,1,4,5,8,9\},$ $\{0,1,4,8,9\},$ and $\{0,1,5,8,9\}.$ It is noteworthy that a larger solution was found before the shorter ones. (This is not invariably the case, though.)

Computing Costs

Storage isn't a problem, because the tree is constructed dynamically, at any stage of a depth-first search only a node and its ascendants are on the stack, and the depth cannot exceed $n-1.$

Computational time is a problem. When the tree is structured like the example, approximately $O(n^2)$ edges need to be traversed. But in other examples, such as $D=\{0,1,2,\ldots,n\},$ the number of solutions grows exponentially with $n,$ implying at least $O(\exp(n))$ edges are traversed. There is thus an unavoidably terrible worst-case computational cost.

Timing studies with random inputs suggest this problem can be efficiently solved when $D$ has been generated from a set with approximately $\sqrt{n}$ elements, because then there tend to be only a small number of solutions. Once the size of $D$ reaches $1.5\sqrt{n},$ the numbers of solutions explode exponentially and therefore so does the time taken to find them all.

Code

Because I have been using R for most of my work for many years, I created a prototype in R (despite its inherent inefficiencies). The code is below. As an example of its output, executing the commands

d <- c(1,6,7,8,12,13,14,15,19,20,21,26,27,33,34)
traverse(energy.new(d))

produces the following output, normalized to match the start-at-1 normalization of the question:

1 2 8 9 15 16 22 28 35
1 2 8 9 16 22 28 35
1 2 8 15 16 22 28 35
1 2 8 16 22 28 35
1 2 9 15 16 22 28 35
1 2 9 16 22 28 35
1 8 14 20 21 27 28 34 35
1 8 14 20 21 27 34 35
1 8 14 20 21 28 34 35
1 8 14 20 27 28 34 35
1 8 14 20 27 34 35
1 8 14 20 28 34 35

These 12 solutions have from $7$ to $9$ elements each. The shortest ones are the four given in the question.

---

#
# (1) Create an "energy" object out of a putative difference set `x`.
#
energy.new <- function(x) {
  if (length(x) > 0) {
    #
    # Normalize the values.
    #
    x <- sort(unique(abs(x)))
    n <- setdiff(max(x), 0)
    if (x[1] == 0) x <- x[-1]
  } else {
    n <- integer(0)
  }
  Delta <- x
  Levels <- n    # Does not include 0
  Spectrum <- n  # Does not include 0
  Candidates <- Delta[(n - Delta) %in% Delta]
  obj <- list(Delta=Delta,
              Levels=Levels,
              Spectrum=Spectrum,
              Candidates=Candidates,
              Complete=isTRUE(all.equal(sort(Delta), sort(Spectrum))))
  class(obj) <- "energy"
  return(obj)
}
#
# (2) Traverse the tree.  Returns the number of solutions found.
#     Calls `emit` to output the solutions.
# Optional arguments allow some control over the traversal and the diagnostic output.
#
traverse <- function(e, stop.at.first=FALSE, shallow=FALSE, verbose=FALSE,
                     depth=0, offset=0) {
  n.solutions <- 0
  if(e$Complete && isTRUE(stop.at.first)) return(n.solutions)

  for(i in seq_along(e$Candidates)) {
    candidate <- e$Candidates[i]
    if(verbose) cat(rep("\t", depth+1),
                    paste0("(",paste(offset + c(0, sort(e$Levels)), collapse=","),")"),
                    "+", offset + candidate, "\n")
    e.0 <- e
    e.0$Candidates <- e.0$Candidates[-seq_len(i)]
    e.0$Levels <- c(e.0$Levels, candidate)
    e.0$Spectrum <-  unique(c(e.0$Spectrum, candidate,
                              abs(candidate - e.0$Levels[e.0$Levels != candidate])))
    e.0$Candidates <- e.0$Candidates[abs(e.0$Candidates - candidate) %in% e.0$Delta]
    if(!isTRUE(e.0$Complete)) {
      e.0$Complete <- isTRUE(all.equal(sort(e.0$Delta), sort(e.0$Spectrum)))
    }
    if(isTRUE(e.0$Complete)) {
      emit(e.0)
      n.solutions <- n.solutions + 1
    }
    if(!e.0$Complete || !isTRUE(shallow)) {
      n.solutions <- n.solutions + traverse(e.0, stop.at.first, shallow, verbose,
                                            depth=depth+1, offset=offset)
    }
    if(n.solutions > 0 && isTRUE(stop.at.first)) break
  }
  return(n.solutions)
}
#------------------------------------------------------------------------------#
#
# Helper functions.
#
spectrum <- function(levels) { # Compute the difference set, without zero.
  rev(sort(unique(abs(c(outer(c(0,levels), c(0,levels), `-`)))))[-1])
}
check <- function(e) {         # Verify a putative solution.
  isTRUE(all.equal(sort(e$Delta), sort(spectrum(e$Levels))))
}
emit <- function(e, offset=1) {# Output a solution.
  cat(paste(offset + c(0, sort(e$Levels))))
  if (check(e)) {
    cat("\n")
  } else {
    cat(" *** Invalid!\n")
  }
}

Answer 2

I was able to generate the same four solutions using a MIP model (written in Java, solved with CPLEX 12.10). The solution time was about 35 ms. on a decent but not spectacular PC. My model uses $O(n)$ binary variables and $O(n^2)$ variables with domain $[0,1]$, where $n$ is the maximum difference, so it might not scale to problems with "large" differences.