diff --git a/Project.toml b/Project.toml
index 3f105e765e..23c18bd6c2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Manopt"
 uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
 authors = ["Ronny Bergmann <manopt@ronnybergmann.net>"]
-version = "0.4.38"
+version = "0.4.39"
 
 [deps]
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
diff --git a/docs/src/solvers/DouglasRachford.md b/docs/src/solvers/DouglasRachford.md
index 396d1aa41b..98c548b858 100644
--- a/docs/src/solvers/DouglasRachford.md
+++ b/docs/src/solvers/DouglasRachford.md
@@ -5,41 +5,45 @@ manifolds in [BergmannPerschSteidl:2016](@cite).
 
 The aim is to minimize the sum
 
-$F(x) = f(x) + g(x)$
+```math
+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
+``\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}(x) = \exp_{\operatorname{prox}_{λ f}(x)} \bigl( -\log_{\operatorname{prox}_{λ f}(x)} x \bigr)$.
+```math
+\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 ∈ \mathbb N} \alpha_k(1-\alpha_k) =  ∈ fty$
-and $λ > 0$ which might depend on iteration $k$ as well) be given.
+Let ``\alpha_k ∈  [0,1]`` with ``\sum_{k ∈ \mathbb N} \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
+Then the (P)DRA algorithm for initial data ``x_0 ∈ \mathcal H`` as
 
 ## Initialization
 
-Initialize $t_0 = x_0$ and $k=0$
+Initialize ``t_0 = x_0`` and ``k=0``
 
 ## Iteration
 
 Repeat until a convergence criterion is reached
 
-1. Compute $s_k = \operatorname{refl}_{λ f}\operatorname{refl}_{λ g}(t_k)$
-2. Within that operation, store $x_{k+1} = \operatorname{prox}_{λ g}(t_k)$ which is the prox the inner reflection reflects at.
-3. Compute $t_{k+1} = g(\alpha_k; t_k, s_k)$
-4. Set $k = k+1$
+1. Compute ``s_k = \operatorname{refl}_{λ f}\operatorname{refl}_{λ g}(t_k)``
+2. Within that operation, store ``p_{k+1} = \operatorname{prox}_{λ g}(t_k)`` which is the prox the inner reflection reflects at.
+3. 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
+4. Set ``k = k+1``
 
 ## Result
 
-The result is given by the last computed $x_K$.
+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
+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),
+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
diff --git a/src/functions/manifold_functions.jl b/src/functions/manifold_functions.jl
index d080f9c51c..e384273476 100644
--- a/src/functions/manifold_functions.jl
+++ b/src/functions/manifold_functions.jl
@@ -33,7 +33,7 @@ retraction, respectively.
 This can also be done in place of `q`.
 
 ## Keyword arguments
-* `retraction_method` - (`default_retration_metiod(M, typeof(p))`) the retraction to use in the reflection
+* `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
diff --git a/src/solvers/DouglasRachford.jl b/src/solvers/DouglasRachford.jl
index 72faca70c0..088370b721 100644
--- a/src/solvers/DouglasRachford.jl
+++ b/src/solvers/DouglasRachford.jl
@@ -11,8 +11,10 @@ Store all options required for the DouglasRachford algorithm,
 * `α` – relaxation of the step from old to new iterate, i.e.
   ``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 inplace or allocating
+* `retraction_method` – a retraction method
 * `stop` – a [`StoppingCriterion`](@ref)
 * `parallel` – indicate whether we are running a parallel Douglas-Rachford or not.
 
@@ -34,7 +36,16 @@ Generate the options for a Manifold `M` and an initial point `p`, where the foll
 * `parallel` – (`false`) indicate whether we are running a parallel Douglas-Rachford
   or not.
 """
-mutable struct DouglasRachfordState{P,Tλ,Tα,TR,S,E} <: AbstractManoptSolverState
+mutable struct DouglasRachfordState{
+    P,
+    Tλ,
+    Tα,
+    TR,
+    S,
+    E<:AbstractEvaluationType,
+    TM<:AbstractRetractionMethod,
+    ITM<:AbstractInverseRetractionMethod,
+} <: AbstractManoptSolverState
     p::P
     p_tmp::P
     s::P
@@ -43,6 +54,8 @@ mutable struct DouglasRachfordState{P,Tλ,Tα,TR,S,E} <: AbstractManoptSolverSta
     α::Tα
     R::TR
     reflection_evaluation::E
+    retraction_method::TM
+    inverse_retraction_method::ITM
     stop::S
     parallel::Bool
     function DouglasRachfordState(
@@ -54,8 +67,19 @@ mutable struct DouglasRachfordState{P,Tλ,Tα,TR,S,E} <: AbstractManoptSolverSta
         reflection_evaluation::E=AllocatingEvaluation(),
         stopping_criterion::S=StopAfterIteration(300),
         parallel=false,
-    ) where {P,Fλ,Fα,FR,S<:StoppingCriterion,E<:AbstractEvaluationType}
-        return new{P,Fλ,Fα,FR,S,E}(
+        retraction_method::TM=default_retraction_method(M, typeof(p)),
+        inverse_retraction_method::ITM=default_inverse_retraction_method(M, typeof(p)),
+    ) where {
+        P,
+        Fλ,
+        Fα,
+        FR,
+        S<:StoppingCriterion,
+        E<:AbstractEvaluationType,
+        TM<:AbstractRetractionMethod,
+        ITM<:AbstractInverseRetractionMethod,
+    }
+        return new{P,Fλ,Fα,FR,S,E,TM,ITM}(
             p,
             copy(M, p),
             copy(M, p),
@@ -64,6 +88,8 @@ mutable struct DouglasRachfordState{P,Tλ,Tα,TR,S,E} <: AbstractManoptSolverSta
             α,
             R,
             reflection_evaluation,
+            retraction_method,
+            inverse_retraction_method,
             stopping_criterion,
             parallel,
         )
@@ -135,13 +161,18 @@ If you provide a [`ManifoldProximalMapObjective`](@ref) `mpo` instead, the proxi
 * `α` – (`(iter) -> 0.9`) relaxation of the step from old to new iterate, i.e.
   ``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)
+* `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
 * `R` – method employed in the iteration to perform the reflection of `x` at the prox `p`.
-  This uses by default `reflect`](@ref) or `reflect!` depending on `reflection_evaluation` and
+  This uses by default [`reflect`](@ref) or `reflect!` depending on `reflection_evaluation` and
   the retraction and inverse retraction specified by `retraction_method` and `inverse_retraction_method`, respectively.
 * `reflection_evaluation` – ([`AllocatingEvaluation`](@ref) whether `R` works inplace or allocating
-* `retraction_method` - (`default_retration_metiod(M, typeof(p))`) the retraction to use in the reflection (ignored, if you set `R` directly)
-* `stopping_criterion` – ([`StopWhenAny`](@ref)`(`[`StopAfterIteration`](@ref)`(200),`[`StopWhenChangeLess`](@ref)`(10.0^-5))`) a [`StoppingCriterion`](@ref).
+* `retraction_method` - (`default_retraction_metiod(M, typeof(p))`) the retraction to use in
+  - the reflection (ignored, if you set `R` directly)
+  - the relaxation step
+* `stopping_criterion` – ([`StopWhenAny`](@ref)`(`[`StopAfterIteration`](@ref)`(200),`[`StopWhenChangeLess`](@ref)`(10.0^-5))`)
+  a [`StoppingCriterion`](@ref).
 * `parallel` – (`false`) clarify that we are doing a parallel DR, i.e. on a
   `PowerManifold` manifold with two proxes. This can be used to trigger
   parallel Douglas–Rachford if you enter with two proxes. Keep in mind, that a
@@ -296,6 +327,8 @@ function DouglasRachford!(
         α=α,
         R=R,
         reflection_evaluation=reflection_evaluation,
+        retraction_method=retraction_method,
+        inverse_retraction_method=inverse_retraction_method,
         stopping_criterion=stopping_criterion,
         parallel=parallel > 0,
     )
@@ -360,7 +393,13 @@ function step_solver!(amp::AbstractManoptProblem, drs::DouglasRachfordState, i)
     get_proximal_map!(amp, drs.p, drs.λ(i), drs.s_tmp, 2)
     _reflect!(M, drs.s_tmp, drs.p, drs.s_tmp, drs.R, drs.reflection_evaluation)
     # relaxation
-    drs.s = shortest_geodesic(M, drs.s, drs.s_tmp, drs.α(i))
+    drs.s = retract(
+        M,
+        drs.s,
+        inverse_retract(M, drs.s, drs.s_tmp, drs.inverse_retraction_method),
+        drs.α(i),
+        drs.retraction_method,
+    )
     return drs
 end
 get_solver_result(drs::DouglasRachfordState) = drs.parallel ? drs.p[1] : drs.p
diff --git a/src/solvers/subgradient.jl b/src/solvers/subgradient.jl
index 473c2326c4..002bf144dd 100644
--- a/src/solvers/subgradient.jl
+++ b/src/solvers/subgradient.jl
@@ -5,7 +5,7 @@ stores option values for a [`subgradient_method`](@ref) solver
 
 # Fields
 
-* `retraction_method` – the retration to use within
+* `retraction_method` – the rectration to use within
 * `stepsize` – ([`ConstantStepsize`](@ref)`(M)`) a [`Stepsize`](@ref)
 * `stop` – ([`StopAfterIteration`](@ref)`(5000)``)a [`StoppingCriterion`](@ref)
 * `p` – (initial or current) value the algorithm is at