Difference-map Algorithm - Algorithm

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 + β ,

f_A(x) = P_A(x) − (P_A(x)−x)/β ,

f_B(x) = P_B(x) + (P_B(x)−x)/β .

The real parameter β 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 β = 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) .

A point x is a fixed point of the map x ↦ D(x) precisely when P_A( f_B(x))= P_B( f_A(x)). 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:

Δ = | P_A( f_B(x)) − P_B( f_A(x)) | .

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