Difference map algorithm

The difference map algorithm is a search algorithm for general constraint satisfaction problems. It is a meta-algorithm in the sense that it is built from more basic algorithms that perform projections onto constraint sets. From a mathematical perspective, the difference map algorithm is a dynamical system based on a mapping of Euclidean space. Solutions are encoded as fixed points of the mapping.

Although originally conceived as a general method for solving the phase problem, the difference map algorithm has been used for the boolean satisfiability problem, protein structure prediction, Ramsey numbers, diophantine equations, and "Sudoku" [V. Elser, I. Rankenburg, and P. Thibault, "Searching with iterated maps". "Proceedings of the National Academy of Sciences" USA. (2007). 104:418-423. http://www.pnas.org/cgi/content/short/104/2/418] . Since these applications include NP-complete problems, the scope of the difference map is that of an incomplete algorithm. Whereas incomplete algorithms can efficiently verify solutions (once a candidate is found), they cannot prove that a solution does not exist.

The difference map algorithm is a generalization of two iterative methods: Fienup's hybrid input-output phase retrieval algorithm [J.R. Fienup, "Phase retrieval algorithms: a comparison". "Applied Optics". (1982). 21:2758-2769.] and the Douglas-Rachford algorithm [H.H. Bauschke, P.L. Combettes, and D.R. Luke, "Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization". "Journal of the Optical Society of America A". (2002). 19:1334-1345.] for convex optimization. Iterative methods, in general, have a long history in phase retrieval and convex optimization. The use of this style of algorithm for hard, non-convex problems is a more recent development.

Algorithm

The problem to be solved must first be formulated as a set intersection problem in Euclidean space: find an "x" in the intersection of sets "A" and "B". Another prerequisite is an implementation of the projections "P"_"A" and "P"_"B" that, given an arbitrary input point "x", return a point in the constraint set "A" or "B" that is nearest to "x". One iteration of the algorithm is given by the mapping:

:"x" → "D"("x") = "x" + β [ "P"_"A"( "f"_"B"("x")) - "P"_"B"( "f"_"A"("x")) ] ,

::"f"_"A"("x") = "P"_"A"("x") - ("P"_"A"("x")-"x")/β ,

::"f"_"B"("x") = "P"_"B"("x") + ("P"_"B"("x")-"x")/β .

The real parameter β can have either sign; optimal values depend on the application and are determined through experimentation. As a first guess, the choice β = 1 (or β = -1) is recommended because it reduces the number of projection computations per iteration:

:"D"("x") = "x" + "P"_"A"(2 "P"_"B"("x") - "x") - "P"_"B"("x") .

The progress of the algorithm is monitored by inspecting the norm of the difference of the two projections:

:"Δ" = | "P"_"A"( "f"_"B"("x")) - "P"_"B"( "f"_"A"("x")) | .

When this vanishes, at fixed points of the map, a point common to both constraint sets has been found and the algorithm is terminated. The set of fixed points in a particular application will normally have a large dimension, even when the solution set is a single point.

Example: logical satisfiability

Incomplete algorithms, such as stochastic local search, are widely used for finding satisfying truth assignments to boolean formulas. As an example of solving an instance of 2-SAT with the difference map algorithm, consider the following formula (~ indicates NOT)::("q"₁ or "q"₂) and (~"q"₁ or "q"₃) and (~"q"₂ or ~"q"₃) and ("q"₁ or ~"q"₂)

To each of the eight literals in this formula we assign one real variable in an eight dimensional Euclidean space. The structure of the 2-SAT formula can be recovered when these variables are arranged in a table:

:

This is followed by 2"P"_B("x"₀) - "x"₀,

:

Here is the second iteration, "D"("x"₁) = "x"₂ :

:

This is a fixed point: "D"("x"₂) = "x"₂. The iterate is unchanged because the two projections agree. From "P"_"B"("x"₂) ,

:

we can read off the satisfying truth assignment: "q"₁ = "TRUE", "q"₂ = "FALSE", "q"₃ = "TRUE".

Chaotic dynamics

In the simple 2-SAT example above, the norm of the difference map increment "Δ" decreased monotonically to zero in three iterations. This contrasts the behavior of "Δ" when the difference map is given a hard instance of 3-SAT, where it fluctuates strongly prior to the discovery of the fixed point. As a dynamical system the difference map is believed to be chaotic, and that the space being searched is a strange attractor.

Phase retrieval

In phase retrieval a signal or image is reconstructed from the modulus (absolute value, magnitude) of its discrete Fourier transform. For example, the source of the modulus data may be the Fraunhofer diffraction pattern formed when an object is illuminated with coherent light.

The projection to the Fourier modulus constraint, say "P"_"A", is accomplished by first computing the discrete Fourier transform of the signal or image, rescaling the moduli to agree with the data, and then inverse transforming the result. This is a projection, in the sense that the Euclidean distance to the constraint is minimized, because (i) the discrete Fourier transform, as a unitary transformation, preserves distance, and (ii) rescaling the modulus (without modifying the phase) is the smallest change that realizes the modulus constraint.

To recover the unknown phases of the Fourier transform the difference map relies on the projection to another constraint, "P"_"B". This may take several forms, as the object being reconstructed may be known to be positive, have a bounded support, etc. In the reconstruction of the Wikipedia logo, for example, the effect of the projection "P"_"B" was to zero all values outside a rectangular support and also to zero all negative values within the support.

References