# Douglas—Rachford algorithm

The (Parallel) Douglas—Rachford ((P)DR) algorithm was generalized to Hadamard manifolds in [BPS16].

The aim is to minimize the sum

\[F(p) = f(p) + g(p)\]

on a manifold, where the two summands have proximal maps $\operatorname{prox}_{λ f}, \operatorname{prox}_{λ g}$ that are easy to evaluate (maybe in closed form, or not too costly to approximate). Further, define the reflection operator at the proximal map as

\[\operatorname{refl}_{λ f}(p) = \operatorname{retr}_{\operatorname{prox}_{λ f}(p)} \bigl( -\operatorname{retr}^{-1}_{\operatorname{prox}_{λ f}(p)} p \bigr).\]

Let $\alpha_k ∈ [0,1]$ with $\sum_{k ∈ ℕ} \alpha_k(1-\alpha_k) = \infty$ and $λ > 0$ (which might depend on iteration $k$ as well) be given.

Then the (P)DRA algorithm for initial data $x_0 ∈ \mathcal H$ as

## Initialization

Initialize $t_0 = x_0$ and $k=0$

## Iteration

Repeat until a convergence criterion is reached

- Compute $s_k = \operatorname{refl}_{λ f}\operatorname{refl}_{λ g}(t_k)$
- Within that operation, store $p_{k+1} = \operatorname{prox}_{λ g}(t_k)$ which is the prox the inner reflection reflects at.
- Compute $t_{k+1} = g(\alpha_k; t_k, s_k)$, where $g$ is a curve approximating the shortest geodesic, provided by a retraction and its inverse
- Set $k = k+1$

## Result

The result is given by the last computed $p_K$.

For the parallel version, the first proximal map is a vectorial version where in each component one prox is applied to the corresponding copy of $t_k$ and the second proximal map corresponds to the indicator function of the set, where all copies are equal (in $\mathcal H^n$, where $n$ is the number of copies), leading to the second prox being the Riemannian mean.

## Interface

`Manopt.DouglasRachford`

— Function```
DouglasRachford(M, f, proxes_f, p)
DouglasRachford(M, mpo, p)
```

Compute the Douglas-Rachford algorithm on the manifold $\mathcal M$, initial data $p$ and the (two) proximal maps `proxMaps`

, see [BPS16].

For $k>2$ proximal maps, the problem is reformulated using the parallel Douglas Rachford: a vectorial proximal map on the power manifold $\mathcal M^k$ is introduced as the first proximal map and the second proximal map of the is set to the `mean`

(Riemannian Center of mass). This hence also boils down to two proximal maps, though each evaluates proximal maps in parallel, that is, component wise in a vector.

If you provide a `ManifoldProximalMapObjective`

`mpo`

instead, the proximal maps are kept unchanged.

**Input**

`M`

: a Riemannian Manifold $\mathcal M$`F`

: a cost function consisting of a sum of cost functions`proxes_f`

: functions of the form`(M, λ, p)-> q`

performing a proximal maps, where`λ`

denotes the proximal parameter, for each of the summands of`F`

. These can also be given in the`InplaceEvaluation`

variants`(M, q, λ p) -> q`

computing in place of`q`

.`p`

: initial data $p ∈ \mathcal M$

**Optional values**

`evaluation`

: (`AllocatingEvaluation`

) specify whether the proximal maps work by allocation (default) form`prox(M, λ, x)`

or`InplaceEvaluation`

in-place`λ`

: (`(iter) -> 1.0`

) function to provide the value for the proximal parameter during the calls`α`

: (`(iter) -> 0.9`

) relaxation of the step from old to new iterate, to be precise $t_{k+1} = g(α_k; t_k, s_k)$, where $s_k$ is the result of the double reflection involved in the DR algorithm`inverse_retraction_method`

- (`default_inverse_retraction_method(M, typeof(p))`

) the inverse retraction to use within- the reflection (ignored, if you set
`R`

directly) - the relaxation step

- the reflection (ignored, if you set
`R`

: method employed in the iteration to perform the reflection of`x`

at the prox`p`

. This uses by default`reflect`

or`reflect!`

depending on`reflection_evaluation`

and the retraction and inverse retraction specified by`retraction_method`

and`inverse_retraction_method`

, respectively.`reflection_evaluation`

: (`AllocatingEvaluation`

whether`R`

works in-place or allocating`retraction_method`

: (`default_retraction_metiod(M, typeof(p))`

) the retraction to use in- the reflection (ignored, if you set
`R`

directly) - the relaxation step

- the reflection (ignored, if you set
`stopping_criterion`

: (`StopAfterIteration`

`(200) |`

`StopWhenChangeLess`

`(1e-5)`

) a`StoppingCriterion`

.`parallel`

: (`false`

) indicate whether to use a parallel Douglas-Rachford or not.

and the ones that are passed to `decorate_state!`

for decorators.

**Output**

the obtained (approximate) minimizer $p^*$, see `get_solver_return`

for details

`DouglasRachford(M, f, proxes_f, p; kwargs...)`

a doc string with some math $t_{k+1} = g(α_k; t_k, s_k)$

`Manopt.DouglasRachford!`

— Function```
DouglasRachford!(M, f, proxes_f, p)
DouglasRachford!(M, mpo, p)
```

Compute the Douglas-Rachford algorithm on the manifold $\mathcal M$, initial data $p ∈ \mathcal M$ and the (two) proximal maps `proxes_f`

in place of `p`

.

For $k>2$ proximal maps, the problem is reformulated using the parallel Douglas Rachford: a vectorial proximal map on the power manifold $\mathcal M^k$ is introduced as the first proximal map and the second proximal map of the is set to the `mean`

(Riemannian Center of mass). This hence also boils down to two proximal maps, though each evaluates proximal maps in parallel, that is component wise in a vector.

While creating the new staring point `p'`

on the power manifold, a copy of `p`

Is created, so that the (by k>2 implicitly generated) parallel Douglas Rachford does not work in-place for now.

If you provide a `ManifoldProximalMapObjective`

`mpo`

instead, the proximal maps are kept unchanged.

**Input**

`M`

: a Riemannian Manifold $\mathcal M$`f`

: a cost function consisting of a sum of cost functions`proxes_f`

: functions of the form`(M, λ, p)->q`

or`(M, q, λ, p)->q`

performing a proximal map, where`λ`

denotes the proximal parameter, for each of the summands of`f`

.`p`

: initial point $p ∈ \mathcal M$

For more options, see `DouglasRachford`

.

## State

`Manopt.DouglasRachfordState`

— Type`DouglasRachfordState <: AbstractManoptSolverState`

Store all options required for the DouglasRachford algorithm,

**Fields**

`p`

: the current iterate (result) For the parallel Douglas-Rachford, this is not a value from the`PowerManifold`

manifold but the mean.`s`

: the last result of the double reflection at the proximal maps relaxed by`α`

.`λ`

: function to provide the value for the proximal parameter during the calls`α`

: relaxation of the step from old to new iterate, to be precise $x^{(k+1)} = g(α(k); x^{(k)}, t^{(k)})$, where $t^{(k)}$ is the result of the double reflection involved in the DR algorithm`inverse_retraction_method`

: an inverse retraction method`R`

: method employed in the iteration to perform the reflection of`x`

at the prox`p`

.`reflection_evaluation`

: whether`R`

works in-place or allocating`retraction_method`

: a retraction method`stop`

: a`StoppingCriterion`

`parallel`

: indicate whether to use a parallel Douglas-Rachford or not.

**Constructor**

`DouglasRachfordState(M, p; kwargs...)`

Generate the options for a Manifold `M`

and an initial point `p`

, where the following keyword arguments can be used

`λ`

: (`(iter)->1.0`

) function to provide the value for the proximal parameter during the calls`α`

: (`(iter)->0.9`

) relaxation of the step from old to new iterate, to be precise $x^{(k+1)} = g(α(k); x^{(k)}, t^{(k)})$, where $t^{(k)}$ is the result of the double reflection involved in the DR algorithm`R`

: (`reflect`

or`reflect!`

) method employed in the iteration to perform the reflection of`x`

at the prox`p`

, which function is used depends on`reflection_evaluation`

.`reflection_evaluation`

: (`AllocatingEvaluation`

`()`

) specify whether the reflection works in-place or allocating (default)`stopping_criterion`

: (`StopAfterIteration`

`(300)`

) a`StoppingCriterion`

`parallel`

: (`false`

) indicate whether to use a parallel Douglas-Rachford or not.

For specific `DebugAction`

s and `RecordAction`

s see also Cyclic Proximal Point.

Furthermore, this solver has a short hand notation for the involved `reflect`

ion.

`Manopt.reflect`

— Function```
reflect(M, f, x; kwargs...)
reflect!(M, q, f, x; kwargs...)
```

reflect the point `x`

from the manifold `M`

at the point `f(x)`

of the function $f: \mathcal M → \mathcal M$, given by

\[ \operatorname{refl}_f(x) = \operatorname{refl}_{f(x)}(x),\]

Compute the result in `q`

.

see also `reflect`

`(M,p,x)`

, to which the keywords are also passed to.

```
reflect(M, p, x, kwargs...)
reflect!(M, q, p, x, kwargs...)
```

Reflect the point `x`

from the manifold `M`

at point `p`

, given by

\[ \operatorname{refl}_p(x) = \operatorname{retr}_p(-\operatorname{retr}^{-1}_p x).\]

where $\operatorname{retr}$ and $\operatorname{retr}^{-1}$ denote a retraction and an inverse retraction, respectively. This can also be done in place of `q`

.

**Keyword arguments**

`retraction_method`

: (`default_retraction_metiod(M, typeof(p))`

) the retraction to use in the reflection`inverse_retraction_method`

: (`default_inverse_retraction_method(M, typeof(p))`

) the inverse retraction to use within the reflection

and for the `reflect!`

additionally

`X`

: (`zero_vector(M,p)`

) a temporary memory to compute the inverse retraction in place. otherwise this is the memory that would be allocated anyways.

```
reflect(M, f, x; kwargs...)
reflect!(M, q, f, x; kwargs...)
```

reflect the point `x`

from the manifold `M`

at the point `f(x)`

of the function $f: \mathcal M → \mathcal M$, given by

\[ \operatorname{refl}_f(x) = \operatorname{refl}_{f(x)}(x),\]

Compute the result in `q`

.

see also `reflect`

`(M,p,x)`

, to which the keywords are also passed to.

```
reflect(M, p, x, kwargs...)
reflect!(M, q, p, x, kwargs...)
```

Reflect the point `x`

from the manifold `M`

at point `p`

, given by

\[ \operatorname{refl}_p(x) = \operatorname{retr}_p(-\operatorname{retr}^{-1}_p x).\]

where $\operatorname{retr}$ and $\operatorname{retr}^{-1}$ denote a retraction and an inverse retraction, respectively. This can also be done in place of `q`

.

**Keyword arguments**

`retraction_method`

: (`default_retraction_metiod(M, typeof(p))`

) the retraction to use in the reflection`inverse_retraction_method`

: (`default_inverse_retraction_method(M, typeof(p))`

) the inverse retraction to use within the reflection

and for the `reflect!`

additionally

`X`

: (`zero_vector(M,p)`

) a temporary memory to compute the inverse retraction in place. otherwise this is the memory that would be allocated anyways.

## Technical details

The `DouglasRachford`

solver requires the following functions of a manifold to be available

- A
`retract!`

`(M, q, p, X)`

; it is recommended to set the`default_retraction_method`

to a favourite retraction. If this default is set, a`retraction_method=`

does not have to be specified. - An
`inverse_retract!`

`(M, X, p, q)`

; it is recommended to set the`default_inverse_retraction_method`

to a favourite retraction. If this default is set, a`inverse_retraction_method=`

does not have to be specified. - A `copyto!
`(M, q, p)`

and`copy`

`(M,p)`

for points.

By default, one of the stopping criteria is `StopWhenChangeLess`

, which requires

- An
`inverse_retract!`

`(M, X, p, q)`

; it is recommended to set the`default_inverse_retraction_method`

to a favourite retraction. If this default is set, a`inverse_retraction_method=`

or`inverse_retraction_method_dual=`

(for $\mathcal N$) does not have to be specified or the`distance`

`(M, p, q)`

for said default inverse retraction.