Difference-map algorithm

Last updated May 06, 2022

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 ,^[1] as well as sphere- and disk-packing problems.^[2] 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 ^[3] and the Douglas-Rachford algorithm ^[4] 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:

{\begin{aligned}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{aligned}}

The real parameter $\beta$ 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 $\beta =1$ (or $\beta =-1$ ) is recommended because it reduces the number of projection computations per iteration:

D(x)=x+P_{A}\left(2P_{B}(x)-x\right)-P_{B}(x)

A point $x$ is a fixed point of the map $x\mapsto D(x)$ precisely when $P_{A}\left(f_{B}(x)\right)=P_{B}\left(f_{A}(x)\right)$ . Since the left-hand side is an element of $A$ and the RHS is an element of $B$ , the equality implies that we have found a common element to the two constraint sets. Note that the fixed point $x$ itself need not belong to either $A$ or $B$ . 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:

\Delta =\left|P_{A}\left(f_{B}(x)\right)-P_{B}\left(f_{A}(x)\right)\right|

.

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, 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:

x₁₁	x₁₂
(x₂₁)		x₂₂
	(x₃₁)	(x₃₂)
x₄₁	(x₄₂)

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₁₁, x₂₁ and x₄₁ correspond to the same boolean variable (q₁) 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₁₁ = -x₂₁ = x₄₁

x₁₂ = -x₃₁ = -x₄₂

x₂₂ = -x₃₂

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₁ = (x₁₁ - x₂₁ + x₄₁) / 3

x₁₁→a₁ x₂₁→ -a₁ x₄₁→a₁

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 3⁴ = 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₀, 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 P_B(x₀) :

1	-1
(1)		-1
	(1)	(-1)
1	(1)

This is followed by 2P_B(x₀) - x₀,

2.5	-1.2
(2.4)		-1.4
	(1.7)	(-1.2)
1.5	(1.9)

and then projected onto the other constraint, P_A(2P_B(x₀) - x₀) :

0.53333	-1.6
(-0.53333)		-0.1
	(1.6)	(0.1)
0.53333	(1.6)

Incrementing x₀ by the difference of the two projections gives the first iteration of the difference map, D(x₀) = x₁ :

-0.96666	-1.4
(-1.93333)		0.3
	(0.9)	(0.3)
0.03333	(0.7)

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

-0.3	-1.4
(-2.6)		-0.7
	(0.9)	(-0.7)
0.7	(0.7)

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

1	-1
(-1)		1
	(1)	(-1)
1	(1)

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 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 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 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

Sudoku Solver - A Sudoku solver based on Difference Map algorithm.

Notes

↑ Elser, V.; Rankenburg, I.; Thibault, P. (9 January 2007). "Searching with iterated maps". Proceedings of the National Academy of Sciences. 104 (2): 418–423. doi: 10.1073/pnas.0606359104 . PMC 1766399 . PMID 17202267.
↑ Gravel, Simon; Elser, Veit (22 September 2008). "Divide and concur: A general approach to constraint satisfaction". Physical Review E. 78 (3): 036706. arXiv: 0801.0222 . Bibcode:2008PhRvE..78c6706G. doi:10.1103/PhysRevE.78.036706. PMID 18851188. S2CID 27814394.
↑ Fienup, J. R. (1 August 1982). "Phase retrieval algorithms: a comparison". Applied Optics. 21 (15): 2758–2769. Bibcode:1982ApOpt..21.2758F. doi:10.1364/AO.21.002758. PMID 20396114. S2CID 10777701 .
↑ Bauschke, Heinz H.; Combettes, Patrick L.; Luke, D. Russell (1 July 2002). "Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization". Journal of the Optical Society of America A. 19 (7): 1334–1345. Bibcode:2002JOSAA..19.1334B. CiteSeerX 10.1.1.75.1070 . doi:10.1364/JOSAA.19.001334. PMID 12095200.

Related Research Articles

In mathematics, the Chinese remainder theorem states that if one knows the remainders of the Euclidean division of an integer n by several integers, then one can determine uniquely the remainder of the division of n by the product of these integers, under the condition that the divisors are pairwise coprime.

In mathematics, the Euclidean algorithm, or Euclid's algorithm, is an efficient method for computing the greatest common divisor (GCD) of two integers (numbers), the largest number that divides them both without a remainder. It is named after the ancient Greek mathematician Euclid, who first described it in his Elements . It is an example of an algorithm, a step-by-step procedure for performing a calculation according to well-defined rules, and is one of the oldest algorithms in common use. It can be used to reduce fractions to their simplest form, and is a part of many other number-theoretic and cryptographic calculations.

Fourier transform Mathematical transform that expresses a function of time as a function of frequency

A Fourier transform (FT) is a mathematical transform that decomposes functions depending on space or time into functions depending on spatial frequency or temporal frequency. An example application would be decomposing the waveform of a musical chord into terms of the intensity of its constituent pitches. The term Fourier transform refers to both the frequency domain representation and the mathematical operation that associates the frequency domain representation to a function of space or time.

The method of least squares is a standard approach in regression analysis to approximate the solution of overdetermined systems by minimizing the sum of the squares of the residuals made in the results of each individual equation.

In mathematics, the Radon transform is the integral transform which takes a function f defined on the plane to a function Rf defined on the (two-dimensional) space of lines in the plane, whose value at a particular line is equal to the line integral of the function over that line. The transform was introduced in 1917 by Johann Radon, who also provided a formula for the inverse transform. Radon further included formulas for the transform in three dimensions, in which the integral is taken over planes. It was later generalized to higher-dimensional Euclidean spaces, and more broadly in the context of integral geometry. The complex analogue of the Radon transform is known as the Penrose transform. The Radon transform is widely applicable to tomography, the creation of an image from the projection data associated with cross-sectional scans of an object.

In mathematics, the Eisenstein integers, occasionally also known as Eulerian integers, are the complex numbers of the form

The Lanczos algorithm is an iterative method devised by Cornelius Lanczos that is an adaptation of power methods to find the $"most useful" eigenvalues and eigenvectors of an Hermitian matrix, where is often but not necessarily much smaller than . Although computationally efficient in principle, the method as initially formulated was not useful, due to its numerical instability.$

In quantum mechanics, the Wigner–Weyl transform or Weyl–Wigner transform is the invertible mapping between functions in the quantum phase space formulation and Hilbert space operators in the Schrödinger picture.

In mathematical optimization, the ellipsoid method is an iterative method for minimizing convex functions. When specialized to solving feasible linear optimization problems with rational data, the ellipsoid method is an algorithm which finds an optimal solution in a number of steps that is polynomial in the input size.

Limited-memory BFGS is an optimization algorithm in the family of quasi-Newton methods that approximates the Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS) using a limited amount of computer memory. It is a popular algorithm for parameter estimation in machine learning. The algorithm's target problem is to minimize $over unconstrained values of the real-vector where is a differentiable scalar function.$

In computational chemistry, a constraint algorithm is a method for satisfying the Newtonian motion of a rigid body which consists of mass points. A restraint algorithm is used to ensure that the distance between mass points is maintained. The general steps involved are: (i) choose novel unconstrained coordinates, (ii) introduce explicit constraint forces, (iii) minimize constraint forces implicitly by the technique of Lagrange multipliers or projection methods.

In algebra, the greatest common divisor of two polynomials is a polynomial, of the highest possible degree, that is a factor of both the two original polynomials. This concept is analogous to the greatest common divisor of two integers.

Phase retrieval is the process of algorithmically finding solutions to the phase problem. Given a complex signal $, of amplitude, and phase :$

Compressed sensing is a signal processing technique for efficiently acquiring and reconstructing a signal, by finding solutions to underdetermined linear systems. This is based on the principle that, through optimization, the sparsity of a signal can be exploited to recover it from far fewer samples than required by the Nyquist–Shannon sampling theorem. There are two conditions under which recovery is possible. The first one is sparsity, which requires the signal to be sparse in some domain. The second one is incoherence, which is applied through the isometric property, which is sufficient for sparse signals.

Spinodal decomposition is a mechanism by which a single thermodynamic phase spontaneously separates into two phases. Decomposition occurs when there is no thermodynamic barrier to phase separation. As a result, phase separation via decomposition does not require the nucleation events resulting from thermodynamic fluctuations which normally trigger phase separation.

The Dirac bracket is a generalization of the Poisson bracket developed by Paul Dirac to treat classical systems with second class constraints in Hamiltonian mechanics, and to thus allow them to undergo canonical quantization. It is an important part of Dirac's development of Hamiltonian mechanics to elegantly handle more general Lagrangians; specifically, when constraints are at hand, so that the number of apparent variables exceeds that of dynamical ones. More abstractly, the two-form implied from the Dirac bracket is the restriction of the symplectic form to the constraint surface in phase space.

Hybrid input-output (HIO) 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 order to properly invert transform the diffraction pattern the phases must be known. Only the amplitude however, can be measured from the intensity of the diffraction pattern and can thus be known experimentally. This fact together with some kind of support constraint can be used in order to iteratively calculate the phases. The HIO algorithm uses negative feedback in Fourier space in order to progressively force the solution to conform to the Fourier domain constraints (support). Unlike the error reduction algorithm which alternately applies Fourier and object constraints the HIO "skips" the object domain step and replaces it with negative feedback acting upon the previous solution.

The multiplicative weights update method is an algorithmic technique most commonly used for decision making and prediction, and also widely deployed in game theory and algorithm design. The simplest use case is the problem of prediction from expert advice, in which a decision maker needs to iteratively decide on an expert whose advice to follow. The method assigns initial weights to the experts, and updates these weights multiplicatively and iteratively according to the feedback of how well an expert performed: reducing it in case of poor performance, and increasing it otherwise. It was discovered repeatedly in very diverse fields such as machine learning, optimization, theoretical computer science, and game theory.

Quantum optimization algorithms are quantum algorithms that are used to solve optimization problems. Mathematical optimization deals with finding the best solution to a problem from a set of possible solutions. Mostly, the optimization problem is formulated as a minimization problem, where one tries to minimize an error which depends on the solution: the optimal solution has the minimal error. Different optimization techniques are applied in various fields such as mechanics, economics and engineering, and as the complexity and amount of data involved rise, more efficient ways of solving optimization problems are needed. The power of quantum computing may allow problems which are not practically feasible on classical computers to be solved, or suggest a considerable speed up with respect to the best known classical algorithm.

Quantum counting algorithm is a quantum algorithm for efficiently counting the number of solutions for a given search problem. The algorithm is based on the quantum phase estimation algorithm and on Grover's search algorithm.

This page is based on this Wikipedia article
Text is available under the CC BY-SA 4.0 license; additional terms may apply.
Images, videos and audio are available under their respective licenses.

[1] Elser, V.; Rankenburg, I.; Thibault, P. (9 January 2007). "Searching with iterated maps". Proceedings of the National Academy of Sciences. 104 (2): 418–423. doi: 10.1073/pnas.0606359104 . PMC 1766399 . PMID 17202267.

[2] Gravel, Simon; Elser, Veit (22 September 2008). "Divide and concur: A general approach to constraint satisfaction". Physical Review E. 78 (3): 036706. arXiv: 0801.0222 . Bibcode:2008PhRvE..78c6706G. doi:10.1103/PhysRevE.78.036706. PMID 18851188. S2CID 27814394.

[3] Fienup, J. R. (1 August 1982). "Phase retrieval algorithms: a comparison". Applied Optics. 21 (15): 2758–2769. Bibcode:1982ApOpt..21.2758F. doi:10.1364/AO.21.002758. PMID 20396114. S2CID 10777701 .

[4] Bauschke, Heinz H.; Combettes, Patrick L.; Luke, D. Russell (1 July 2002). "Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization". Journal of the Optical Society of America A. 19 (7): 1334–1345. Bibcode:2002JOSAA..19.1334B. CiteSeerX 10.1.1.75.1070 . doi:10.1364/JOSAA.19.001334. PMID 12095200.

[1]

[2]

[3]

[4]