Edit on GitHub

AlternatingProjections.jl

A package implementing Alternating Projections methods in Julia.

Alternating projections (theoretical background)

Alternating Projections (AP) is a simple method of solving a feasibility problem

\[\text{for sets } A, B, \text{ find } x \in A \cap B\]

if $A\cap B \neq \varnothing .$ In case of inconsistent problem $A\cap B = \varnothing$, AP finds $x \in A$ closest to $B$ in some sense.

The algorithm starts with some initial value $x^0$ and proceeds by projecting $x^k$ alternatively on $A$ and $B$ :

\[y^{k} = P_B(x^k)\]

and

\[x^{k+1} = P_A(y^k).\]

The method was originally proposed by von Neumann in 1949 for $A$ and $B$ being linear subspaces, and was later generalised for any number of convex sets $A_1, \ldots,A_N$ as

\[x^{k+1} = P_{A_{n(k)}}(x^k). \]

It was also recently shown that the method can also work for not convex sets, if the condition of transversality is satisfied, (see for instance the Gerchberg-Saxton method).

A lot of popular algorithms can be explained in AP framework, even if the original algorithm was invented based on other principles. See, for instance, Gerchberg-Saxton algorithm for phase retrieval.

Projections

Projection operator is often more easily described in math formula than in a computer program. It is easy to project on a (multidimensional) sphere, for instance, or to a polyhedron. The projection on a general convex set can be much more difficult (consider, for instance, projection on an ellipse). This explains why in this package we limit ourselves to special types of sets that allows relatively easiness of defining a projection operator.

Forward and backward operators

In some case the projection can be more easily computed in a transformed space (e.g. in the Fourier domain). This can be generalised further by introducing some abstract (even not necassarily linear) forward and backward transforms

\[x^{k+1} = b_{n(k)}(P_{A_{n(k)}}(f_{n(k)}(x^k))), \]

where $f_{n(k)}(x)$ and $b_{n(k)}(x)$ are the forward and backward transforms correspondingly. For instance, in TIP algorithm $f(x) = b(x) = i /_* x$, where $/_*$ is ''deconvolution'' operator, such that $x /_* x = \delta_0$.

Package features

  • types for the frequently used feasible sets and projections on them
  • a general AP algorithm and examples of its adaptation to popular AP algorithms, including:[1]
    • Gerchberg-Saxton
    • Vector Gerchberg-Saxton
    • Gerchberg-Papoulis
    • TIP
    • Douglas–Rachford
    • DRAP

Manual outline

Function Index

AlternatingProjections

The types defined here are absract and not related exclusevely to the alternating projections methos.

AlternatingProjections — Module

AlternatingProjections

Module contains types and functions for formulation and solving feasibility problem using projections-based methods.

source
AlternatingProjections.backward! — Method
backward!(ts::TransformedSet)

Two-argument forward transform assosiated with the set `ts, ts = F(s): p ∈ s ⇔ q = F(p) ∈ ts, it updates the first argument: use backward!(ts)(p,q) to obtain p: q = F(p).

source
AlternatingProjections.forward! — Method
forward!(ts::TransformedSet)

Two-argument forward transform assosiated with the set ts, ts = F(s): p ∈ s ⇔ q = F(p) ∈ ts, it updates the first argument: use forward!(ts)(q,p) to obtain q = F(p).

source
AlternatingProjections.ACset — Type
ACset{T,N,M,K}(amp,projdims)

Constructs an extended version of the AmplitudeCosntrained set. Here, amp can be a tuple of amplitudes, and projdims are the dimeshoins, along with the length of the vector is measured.

TODO extend description

source
AlternatingProjections.AP — Type
AP <: ProjectionsMethod

Classical Alternating Projections method

Project on one set, than on another, stop after fixed number of iterations or if the desired accuracy is achieved.

source
AlternatingProjections.ConstrainedByShape — Type

ConstrainedByShape{T,N}

Set of abstract arryas with elemet typ T and dimensions N defined by the shape constraint|x| = s ampfor somes`. Field n contains the sum of square of all elements.

source
AlternatingProjections.ConstrainedByShape — Method

ConstrainedByShape(a::AbstractArray{T,N})

Construct set defined by the shape constraint |x| = s a for some scaling s. Type and dimension of the set are inhereited from the array.

source
AlternatingProjections.ConstrainedByShapeSaturated — Type

ConstrainedByShapeSaturated(a, mask::Vector) ConstrainedByShapeSaturated(a, AbstractArray{Bool})

The amplitude constrained set only in the indexes given by mask: |xᵢ| = s aᵢ for i ∈ mask for some s. Outside the mask function should be larger than the satruation level |xᵢ| > s for i ∉ mask.

For this set, the amplitude should be provided in the range from 0 to 1 (1 corresponding to the saturated values).

source
AlternatingProjections.ConstrainedBySupport — Type
ConstrainedBySupport(support)

Special type of convex set.

For continuous case: consists of all functions that equals zero outside some fixed area called support: $𝓐_S = \{f : f(x) = 0, x ∉ S \} .$

For discrete case: all arrays that equals zero for indexes outside some index set: $𝓐_S = \{x : x[i] = 0, i ∉ S \} .$

Currently supports only discrete case, with the support defined as a boolean array.

Examples

julia> S =  ConstrainedBySupport([true, false,true])
ConstrainedBySupport(Bool[1, 0, 1])

julia> x = [1, 2, 3]; project(x, S)
3-element Vector{Int64}:
 1
 0
 3

julia> S = ConstrainedBySupport([x^2 + y^2 <=1  for x in -2:.5:2, y in -2:.5:2]);

julia> x = ones(size(S.support));

julia> project(x,S)
9×9 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
source
AlternatingProjections.ConstrainedBySupportNormed — Type
ConstrainedBySupportNormed(support, norm)

ConstrainedBySupportNormed(support, norm) is a set represented by intersection of the ConstrainedBySupport(support) set and a hypersphere of radius norm.

Examples

julia>  a = rand((1,10),(5,5));

julia> mask = a .> 5;

julia> aset= ConstrainedBySupportNormed(mask, 10);

julia> b = project(Float64.(a), aset);

julia> sum(abs2,b) ≈ 10^2
true
source
AlternatingProjections.DR — Type
DR <: ProjectionsMethod

Douglas-Rachford method

Reflect with respect one set, than to another and move there halfway: x → (x + R₂R₁x)/2. Stop after fixed number of iterations or if the desired accuracy is achieved

source
AlternatingProjections.DRAP — Type
DRAP <: ProjectionsMethod

Douglas-Rachford and Alternating Projection method

Combination of DR and AP method (see [1]). Stop after fixed number of iterations or if the desired accuracy is achieved

[1] Nguyen Hieu Thao, Oleg Soloviev, and Michel Verhaegen. Convex combination of alternating projection and Douglas-Rachford operators for phase retrieval. Adv. Comput. Math.

source
AlternatingProjections.FeasibleSet — Type
FeasibleSet

Abstract type representing a set for feasibility problem. For any set we should be able to take it's representative element.

For any concrete subtype, we define a project operator (maybe set-valued).

source
AlternatingProjections.Opiterator1 — Type
Opiterator1(x⁰, f)

Create a simple iteration with operator f and inital value x⁰.

Examples

julia> ttt = Opiterator([3], x -> x .+ 2)
AlternatingProjections.Opiterator1{Vector{Int64}}([3], var"#11#12"())

julia> for state in Base.Iterators.take(ttt,5)
       println(state.xᵏ[1])
       end
5
7
9
11
13

julia> ttt20 = Iterators.take(ttt,20);

julia> for state in ttt20
       print(state.xᵏ[1])
       end
5791113151719212325272931333537394143

Example with Images.jl

julia> using Images

julia> img = Gray.(rand(Float64, (256,256)))

julia> blurop = Opiterator(img, x -> imfilter(x, Kernel.gaussian(2)));

julia> blurop20 = Iterators.take(blurop, 20);

julia> mosaic([copy(b.xᵏ⁻¹) for b in blurop20]; nrow=4, rowmajor = true)
source
AlternatingProjections.ProjectionsMethod — Type

ProjectionsMethod is class of iterative algorithms for solving feasibility problems based on projections on the feasible sets. Examples are Alternating Projections (AP), also known in the literature as Projections on the Convex Sets (POCS), Douglas-Rachford algorithm (DR), their combination DRAP and many others.

source

AbstractProblems

The types defined here are absract and not related exclusevely to the alternating projections methos.

AlternatingProjections.AbstractProblems.solve — Method
solve(p::Problem,x⁰,alg::Algorithm)

solve(p::Problem,alg::IterativeAlgorithm; x⁰, ϵ, maxit, keephistory, snapshshots)

solve(p::Problem,(alg1, alg2,...); x⁰, ϵ, maxit, keephistory, snapshshots)

Solve problem p, using method alg. For iterative algorithms the arguments may be specified separately. Optionally keep the error history and the iteration snapshots.

If sequence of the iterative algorithms is given, they are run subsequently, using the last state of the previous algorith as the inital value.

source
AlternatingProjections.AbstractProblems.IterativeAlgorithm — Type

Iterative methods form an important class of algorithms with iteratively adjust solution.

Concrete types of the IterativeAlgorithm should contain the initial value, tolerance and maximum number of iterations and are used more for convenience. These can be obtained by fucntions initial, tolerance, maxit fuctions. If applied the abstract types, these fucntions return missing and can trigger using of the default values.

In addition, the concrete types contain instructions on whether or not to keeep the history of the convergence and some snapshots of the inner state of an iterator. These values are given by functions keephistory and snapshots.

source

Index

  • 1not implemented features are shown with italics