Difference map algorithm
Encyclopedia
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
, as well as sphere- and disk-packing problems. 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
and the Douglas-Rachford algorithm 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.
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:
The progress of the algorithm is monitored by inspecting the norm of the difference of the two projections:
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.
, 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): and (~q1 or q3) and (~q2 or ~q3) and (q1 or ~q2)
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:
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 x11, x21 and x41 correspond to the same boolean variable (q1) 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:
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:
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:
→ (-1, 1) → (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 x0, say
Using β = 1, the next step is to compute PB(x0) :
This is followed by 2PB(x0) - x0,
and then projected onto the other constraint, PA(2PB(x0) - x0) :
Incrementing x0 by the difference of the two projections gives the first iteration of the difference map, D(x0) = x1 :
Here is the second iteration, D(x1) = x2 :
This is a fixed point: D(x2) = x2. The iterate is unchanged because the two projections agree. From PB(x2) ,
we can read off the satisfying truth assignment: q1 = TRUE, q2 = FALSE, q3 = TRUE.
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.
(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 PA, 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, PB. 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 surface image, for example, the effect of the projection PB was to nullify all values outside a rectangular support, and also to nullify all negative values within the support.
Search algorithm
In computer science, a search algorithm is an algorithm for finding an item with specified properties among a collection of items. The items may be stored individually as records in a database; or may be elements of a search space defined by a mathematical formula or procedure, such as the roots...
for general constraint satisfaction
Constraint satisfaction
In artificial intelligence and operations research, constraint satisfaction is the process of finding a solution to a set of constraints that impose conditions that the variables must satisfy. A solution is therefore a vector of variables that satisfies all constraints.The techniques used in...
problems. It is a meta-algorithm in the sense that it is built from more basic algorithms that perform projections
Projection (linear algebra)
In linear algebra and functional analysis, a projection is a linear transformation P from a vector space to itself such that P2 = P. It leaves its image unchanged....
onto constraint
Constraint (mathematics)
In mathematics, a constraint is a condition that a solution to an optimization problem must satisfy. There are two types of constraints: equality constraints and inequality constraints...
sets. From a mathematical perspective, the difference-map algorithm is a dynamical system
Dynamical system
A dynamical system is a concept in mathematics where a fixed rule describes the time dependence of a point in a geometrical space. Examples include the mathematical models that describe the swinging of a clock pendulum, the flow of water in a pipe, and the number of fish each springtime in a...
based on a mapping
Map (mathematics)
In most of mathematics and in some related technical fields, the term mapping, usually shortened to map, is either a synonym for function, or denotes a particular kind of function which is important in that branch, or denotes something conceptually similar to a function.In graph theory, a map is a...
of Euclidean space
Euclidean space
In mathematics, Euclidean space is the Euclidean plane and three-dimensional space of Euclidean geometry, as well as the generalizations of these notions to higher dimensions...
. Solutions are encoded as fixed points
Fixed point (mathematics)
In mathematics, a fixed point of a function is a point that is mapped to itself by the function. A set of fixed points is sometimes called a fixed set...
of the mapping.
Although originally conceived as a general method for solving the phase problem
Phase problem
In physics the phase problem is the name given to the problem of loss of information concerning the phase that can occur when making a physical measurement. The name itself comes from the field of x-ray crystallography, where the phase problem has to be solved for the determination of a structure...
, the difference-map algorithm has been used for the boolean satisfiability problem
Boolean satisfiability problem
In computer science, satisfiability is the problem of determining if the variables of a given Boolean formula can be assigned in such a way as to make the formula evaluate to TRUE...
, protein structure prediction
Protein structure prediction
Protein structure prediction is the prediction of the three-dimensional structure of a protein from its amino acid sequence — that is, the prediction of its secondary, tertiary, and quaternary structure from its primary structure. Structure prediction is fundamentally different from the inverse...
, Ramsey numbers, diophantine equations, and Sudoku
Sudoku
is a logic-based, combinatorial number-placement puzzle. The objective is to fill a 9×9 grid with digits so that each column, each row, and each of the nine 3×3 sub-grids that compose the grid contains all of the digits from 1 to 9...
, as well as sphere- and disk-packing problems. Since these applications include NP-complete
NP-complete
In computational complexity theory, the complexity class NP-complete is a class of decision problems. A decision problem L is NP-complete if it is in the set of NP problems so that any given solution to the decision problem can be verified in polynomial time, and also in the set of NP-hard...
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
Hybrid input output (HIO) algorithm for phase retrieval
Hybrid input-output algorithm for phase retrieval is a modification of the error reduction algorithm for retrieving the phases in Coherent diffraction imaging. Determining the phases of a diffraction pattern is crucial since the diffraction pattern of an object is its Fourier transform and in...
and the Douglas-Rachford algorithm 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 PA and PB 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 + β [ PA( fB(x)) - PB( fA(x)) ] ,
-
- fA(x) = PA(x) - (PA(x)-x)/β ,
-
- fB(x) = PB(x) + (PB(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 + PA(2 PB(x) - x) - PB(x) .
The progress of the algorithm is monitored by inspecting the norm of the difference of the two projections:
- Δ = | PA( fB(x)) - PB( fA(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
Dimension
In physics and mathematics, the dimension of a space or object is informally defined as the minimum number of coordinates needed to specify any point within it. Thus a line has a dimension of one because only one coordinate is needed to specify a point on it...
, even when the solution set is a single point.
Example: logical satisfiability
Incomplete algorithms, such as stochastic local searchLocal search (optimization)
In computer science, local search is a metaheuristic method for solving computationally hard optimization problems. Local search can be used on problems that can be formulated as finding a solution maximizing a criterion among a number of candidate solutions...
, 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): and (~q1 or q3) and (~q2 or ~q3) and (q1 or ~q2)
To each of the eight literals
Literal (mathematical logic)
In mathematical logic, a literal is an atomic formula or its negation.The definition mostly appears in proof theory , e.g...
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:
| x11 | x12 | |
| (x21) | x22 | |
| (x31) | (x32) | |
| x41 | (x42) |
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 x11, x21 and x41 correspond to the same boolean variable (q1) 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:
- x11 = -x21 = x41
- x12 = -x31 = -x42
- x22 = -x32
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:
- a1 = (x11 - x21 + x41) / 3
- x11 → a1 x21 → -a1 x41 → a1
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:
→ (-1, 1) → (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 x0, say
| -0.5 | -0.8 | |
| (-0.4) | -0.6 | |
| (0.3) | (-0.8) | |
| 0.5 | (0.1) |
Using β = 1, the next step is to compute PB(x0) :
| 1 | -1 | |
| (1) | -1 | |
| (1) | (-1) | |
| 1 | (1) |
This is followed by 2PB(x0) - x0,
| 2.5 | -1.2 | |
| (2.4) | -1.4 | |
| (1.7) | (-1.2) | |
| 1.5 | (1.9) |
and then projected onto the other constraint, PA(2PB(x0) - x0) :
| 0.53333 | -1.6 | |
| (-0.53333) | -0.1 | |
| (1.6) | (0.1) | |
| 0.53333 | (1.6) |
Incrementing x0 by the difference of the two projections gives the first iteration of the difference map, D(x0) = x1 :
| -0.96666 | -1.4 | |
| (-1.93333) | 0.3 | |
| (0.9) | (0.3) | |
| 0.03333 | (0.7) |
Here is the second iteration, D(x1) = x2 :
| -0.3 | -1.4 | |
| (-2.6) | -0.7 | |
| (0.9) | (-0.7) | |
| 0.7 | (0.7) |
This is a fixed point: D(x2) = x2. The iterate is unchanged because the two projections agree. From PB(x2) ,
| 1 | -1 | |
| (-1) | 1 | |
| (1) | (-1) | |
| 1 | (1) |
we can read off the satisfying truth assignment: q1 = TRUE, q2 = FALSE, q3 = TRUE.
Chaotic dynamics
Chaos theory
Chaos theory is a field of study in mathematics, with applications in several disciplines including physics, economics, biology, and philosophy. Chaos theory studies the behavior of dynamical systems that are highly sensitive to initial conditions, an effect which is popularly referred to as the...
, and that the space being searched is a strange attractor.
Phase retrieval
In phase retrieval a signal or image is reconstructed from the modulusAbsolute value
In mathematics, the absolute value |a| of a real number a is the numerical value of a without regard to its sign. So, for example, the absolute value of 3 is 3, and the absolute value of -3 is also 3...
(absolute value, magnitude) of its discrete Fourier transform
Discrete Fourier transform
In mathematics, the discrete Fourier transform is a specific kind of discrete transform, used in Fourier analysis. It transforms one function into another, which is called the frequency domain representation, or simply the DFT, of the original function...
. For example, the source of the modulus data may be the Fraunhofer diffraction
Fraunhofer diffraction
In optics, the Fraunhofer diffraction equation is used to model the diffraction of waves when the diffraction pattern is viewed at a long distance from the diffracting object, and also when it is viewed at the focal plane of an imaging lens....
pattern formed when an object is illuminated with coherent light.
The projection to the Fourier modulus constraint, say PA, 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
Unitary transformation
In mathematics, a unitary transformation may be informally defined as a transformation that respects the inner product: the inner product of two vectors before the transformation is equal to their inner product after the 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, PB. This may take several forms, as the object being reconstructed may be known to be positive, have a bounded support
Support (mathematics)
In mathematics, the support of a function is the set of points where the function is not zero, or the closure of that set . This concept is used very widely in mathematical analysis...
, etc. In the reconstruction of the surface image, for example, the effect of the projection PB was to nullify all values outside a rectangular support, and also to nullify all negative values within the support.