imported>Kvng: add link

2025-06-16T16:30:12Z

add link

New page

[[File:Iterations 0, 100, 200, 300 and 400 in difference-map reconstruction of [[grayscale image]] from Fourier transform modulus.png|thumb|right|Iterations 0, 100, 200, 300 and 400 in the difference-map reconstruction of a grayscale image from its Fourier transform modulus]]

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 [[Projection (linear algebra)|projections]] onto [[Constraint (mathematics)|constraint]] sets. From a mathematical perspective, the difference-map algorithm is a [[dynamical system]] based on a [[Map (mathematics)|mapping]] of [[Euclidean space]]. Solutions are encoded as [[Fixed point (mathematics)|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]]'',<ref>{{cite journal |last1=Elser |first1=V. |last2=Rankenburg |first2=I. |last3=Thibault |first3=P. |title=Searching with iterated maps |journal=Proceedings of the National Academy of Sciences |date=9 January 2007 |volume=104 |issue=2 |pages=418–423 |doi=10.1073/pnas.0606359104 |pmid=17202267 |pmc=1766399 |doi-access=free }}</ref> as well as sphere- and disk-packing problems.<ref>{{cite journal |last1=Gravel |first1=Simon |last2=Elser |first2=Veit |title=Divide and concur: A general approach to constraint satisfaction |journal=Physical Review E |date=22 September 2008 |volume=78 |issue=3 |pages=036706 |doi=10.1103/PhysRevE.78.036706 |pmid=18851188 |arxiv=0801.0222 |bibcode=2008PhRvE..78c6706G |s2cid=27814394 }}</ref> 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 (HIO) algorithm for phase retrieval]]<ref>{{cite journal |last1=Fienup |first1=J. R. |title=Phase retrieval algorithms: a comparison |journal=Applied Optics |date=1 August 1982 |volume=21 |issue=15 |pages=2758–2769 |doi=10.1364/AO.21.002758 |pmid=20396114 |bibcode=1982ApOpt..21.2758F |s2cid=10777701 |s2cid-access=free }}</ref> and the [[Douglas-Rachford algorithm]]<ref>{{cite journal |last1=Bauschke |first1=Heinz H. |last2=Combettes |first2=Patrick L. |last3=Luke |first3=D. Russell |title=Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization |journal=Journal of the Optical Society of America A |date=1 July 2002 |volume=19 |issue=7 |pages=1334–1345 |doi=10.1364/JOSAA.19.001334 |pmid=12095200 |bibcode=2002JOSAA..19.1334B |citeseerx=10.1.1.75.1070 }}</ref> 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 <math>x</math> in the intersection of sets <math>A</math> and <math>B</math>. Another prerequisite is an implementation of the projections <math>P_A</math> and <math>P_B</math> that, given an arbitrary input point <math>x</math>, return a point in the constraint set <math>A</math> or <math>B</math> that is nearest to <math>x</math>. One iteration of the algorithm is given by the mapping:

:<math>
\begin{align}
x \mapsto D(x) &= x + \beta \left[ P_A \left( f_B(x)\right) - P_B \left( f_A(x)\right)\right], \\
f_A(x) &= P_A(x) - \frac{1}{\beta}\left( P_A(x) - x\right), \\
f_B(x) &= P_B(x) + \frac{1}{\beta}\left( P_B(x) - x\right)
\end{align}
</math>
The real parameter <math>\beta</math> should not be equal to 0 but can have either sign; optimal values depend on the application and are determined through experimentation. As a first guess, the choice <math>\beta = 1</math> (or <math>\beta = -1</math>) is recommended because it reduces the number of projection computations per iteration:

:<math>D(x) = x + P_A\left( 2P_B(x) - x\right)-P_B(x)</math>

A point <math>x</math> is a fixed point of the map <math>x \mapsto D(x)</math> precisely when <math>P_A\left(f_B(x)\right) = P_B\left(f_A(x)\right)</math>. Since the left-hand side is an element of <math>A</math> and the RHS is an element of <math>B</math>, the equality implies that we have found a common element to the two constraint sets. Note that the fixed point <math>x</math> itself need not belong to either <math>A</math> or <math>B</math>. The set of fixed points will typically have much higher dimension than the set of solutions.

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

:<math>\Delta = \left| P_A \left( f_B(x)\right) - P_B \left( f_A(x)\right)\right|</math>.

When this vanishes, a point common to both constraint sets has been found and the algorithm can be terminated.

==Example: logical satisfiability ==
Incomplete algorithms, such as stochastic [[Local search (optimization)|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''1 or ''q''2) and (~''q''1 or ''q''3) and (~''q''2 or ~''q''3) and (''q''1 or ~''q''2)

To each of the eight [[Literal (mathematical logic)|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:

:{| class="wikitable" style="background-color:white;text-align: center;"
|-
|width="80pt"|''x''11
|width="80pt"|''x''12
|width="80pt"|
|-
|(''x''21)
|
|''x''22
|-
|
|(''x''31)
|(''x''32)
|-
|''x''41
|(''x''42)
|
|}

Rows are the clauses in the 2-SAT formula and literals corresponding to the same boolean variable are arranged in columns, with negation indicated by parentheses. For example, the real variables ''x''11, ''x''21 and ''x''41 correspond to the same boolean variable (''q''1) or its negation, and are called '''replicas'''.
It is convenient to associate the values 1 and -1 with ''TRUE'' and ''FALSE'' rather than the traditional 1 and 0. With this convention, the compatibility between the replicas takes the form of the following linear equations:

:''x''11 = -''x''21 = ''x''41
:''x''12 = -''x''31 = -''x''42
:''x''22 = -''x''32

The linear subspace where these equations are satisfied is one of the constraint spaces, say ''A'', used by the difference map. To project to this constraint we replace each replica by the signed replica average, or its negative:

:''a''1 = (''x''11 - ''x''21 + ''x''41) / 3
:''x''11 → ''a''1   ''x''21 → -''a''1   ''x''41 → ''a''1

The second difference-map constraint applies to the rows of the table, the clauses. In a satisfying assignment, the two variables in each row must be assigned the values (1, 1), (1, -1), or (-1, 1). The corresponding constraint set, ''B'', is thus a set of 34 = 81 points. In projecting to this constraint the following operation is applied to each row. First, the two real values are rounded to 1 or -1; then, if the outcome is (-1, -1), the larger of the two original values is replaced by 1. Examples:

:(-.2, 1.2) → (-1, 1)
:(-.2, -.8) → (1, -1)

It is a straightforward exercise to check that both of the projection operations described minimize the Euclidean distance between input and output values. Moreover, if the algorithm succeeds in finding a point ''x'' that lies in both constraint sets, then we know that (i) the clauses associated with ''x'' are all ''TRUE'', and (ii) the assignments to the replicas are consistent with a truth assignment to the original boolean variables.

To run the algorithm one first generates an initial point ''x''0, say

:{| class="wikitable" style="background-color:white;text-align: center;"
|-
|width="80pt"|-0.5
|width="80pt"|-0.8
|width="80pt"|
|-
|(-0.4)
|
| -0.6
|-
|
|(0.3)
|(-0.8)
|-
|0.5
|(0.1)
|
|}

Using β = 1, the next step is to compute ''P''B(''x''0) :

:{| class="wikitable" style="background-color:white;text-align: center;"
|-
|width="80pt"|1
|width="80pt"|-1
|width="80pt"|
|-
|(1)
|
| -1
|-
|
|(1)
|(-1)
|-
|1
|(1)
|
|}

This is followed by 2''P''B(''x''0) - ''x''0,

:{| class="wikitable" style="background-color:white;text-align: center;"
|-
|width="80pt"|2.5
|width="80pt"|-1.2
|width="80pt"|
|-
|(2.4)
|
| -1.4
|-
|
|(1.7)
|(-1.2)
|-
|1.5
|(1.9)
|
|}

and then projected onto the other constraint, ''P''A(2''P''B(''x''0) - ''x''0) :

:{| class="wikitable" style="background-color:white;text-align: center;"
|-
|width="80pt"|0.53333
|width="80pt"|-1.6
|width="80pt"|
|-
|(-0.53333)
|
| -0.1
|-
|
|(1.6)
|(0.1)
|-
|0.53333
|(1.6)
|
|}

Incrementing ''x''0 by the difference of the two projections gives the first iteration of the difference map, ''D''(''x''0) = ''x''1 :

:{| class="wikitable" style="background-color:white;text-align: center;"
|-
|width="80pt"|-0.96666
|width="80pt"|-1.4
|width="80pt"|
|-
|(-1.93333)
|
| 0.3
|-
|
|(0.9)
|(0.3)
|-
|0.03333
|(0.7)
|
|}

Here is the second iteration, ''D''(''x''1) = ''x''2 :

:{| class="wikitable" style="background-color:white;text-align: center;"
|-
|width="80pt"|-0.3
|width="80pt"|-1.4
|width="80pt"|
|-
|(-2.6)
|
| -0.7
|-
|
|(0.9)
|(-0.7)
|-
|0.7
|(0.7)
|
|}

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

:{| class="wikitable" style="background-color:white;text-align: center;"
|-
|width="80pt"|1
|width="80pt"|-1
|width="80pt"|
|-
|(-1)
|
| 1
|-
|
|(1)
|(-1)
|-
|1
|(1)
|
|}

we can read off the satisfying truth assignment: ''q''1 = ''TRUE'', ''q''2 = ''FALSE'', ''q''3 = ''TRUE''.

== Chaotic dynamics ==

[[Image:Time series of norm of difference-map increment Δ, during solving random 3-SAT instance.png|right|thumb|Time series of the norm of the difference-map increment ''Δ'' in the course of solving a random 3-SAT instance with 1000 variables and 4200 clauses.]]

In the simple 2-SAT example above, the norm of the difference-map increment ''Δ'' decreased monotonically to zero in three iterations. This contrasts with 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 [[Chaos theory|chaotic]], and that the space being searched is a [[strange attractor]].

== Phase retrieval ==

[[Image:Diffraction data.png|left|thumb|Fourier transform modulus (diffraction pattern) of the grayscale image shown being reconstructed at the top of the page.]]

In phase retrieval a signal or image is reconstructed from the [[Absolute value|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 (mathematics)|support]], etc. In the reconstruction of the surface image, for example, the effect of the projection ''P''''B'' was to nullify all values outside a rectangular support, and also to nullify all negative values within the support.
{{-}}

== External links ==
* [https://barisdemiroz.github.io/sudoku/sudokusolver_demo.html Sudoku Solver] - A Sudoku solver based on Difference Map algorithm.

==Notes==

{{reflist|2}}

[[Category:Search algorithms]]
[[Category:Constraint programming]]

Difference-map algorithm - Revision history

imported>Kvng: add link