The Riemannian Chambolle-Pock algorithm
The Riemannian Chambolle—Pock is a generalization of the Chambolle—Pock algorithm Chambolle and Pock [CP11] It is also known as primal-dual hybrid gradient (PDHG) or primal-dual proximal splitting (PDPS) algorithm.
In order to minimize over $p∈\mathcal M$ the cost function consisting of In order to minimize a cost function consisting of
\[F(p) + G(Λ(p)),\]
over $p∈\mathcal M$
where $F:\mathcal M → \overline{ℝ}$, $G:\mathcal N → \overline{ℝ}$, and $Λ:\mathcal M →\mathcal N$. If the manifolds $\mathcal M$ or $\mathcal N$ are not Hadamard, it has to be considered locally only, that is on geodesically convex sets $\mathcal C \subset \mathcal M$ and $\mathcal D \subset\mathcal N$ such that $Λ(\mathcal C) \subset \mathcal D$.
The algorithm is available in four variants: exact versus linearized (see variant
) as well as with primal versus dual relaxation (see relax
). For more details, see Bergmann, Herzog, Silva Louzeiro, Tenbrinck and Vidal-Núñez [BHS+21]. In the following description is the case of the exact, primal relaxed Riemannian Chambolle—Pock algorithm.
Given base points $m∈\mathcal C$, $n=Λ(m)∈\mathcal D$, initial primal and dual values $p^{(0)} ∈\mathcal C$, $ξ_n^{(0)} ∈T_n^*\mathcal N$, and primal and dual step sizes $\sigma_0$, $\tau_0$, relaxation $\theta_0$, as well as acceleration $\gamma$.
As an initialization, perform $\bar p^{(0)} \gets p^{(0)}$.
The algorithms performs the steps $k=1,…,$ (until a StoppingCriterion
is fulfilled with)
\[ξ^{(k+1)}_n = \operatorname{prox}_{\tau_k G_n^*}\Bigl(ξ_n^{(k)} + \tau_k \bigl(\log_n Λ (\bar p^{(k)})\bigr)^\flat\Bigr)\]
\[p^{(k+1)} = \operatorname{prox}_{\sigma_k F}\biggl(\exp_{p^{(k)}}\Bigl( \operatorname{PT}_{p^{(k)}\gets m}\bigl(-\sigma_k DΛ(m)^*[ξ_n^{(k+1)}]\bigr)^\sharp\Bigr)\biggr)\]
- Update
- $\theta_k = (1+2\gamma\sigma_k)^{-\frac{1}{2}}$
- $\sigma_{k+1} = \sigma_k\theta_k$
- $\tau_{k+1} = \frac{\tau_k}{\theta_k}$
\[\bar p^{(k+1)} = \exp_{p^{(k+1)}}\bigl(-\theta_k \log_{p^{(k+1)}} p^{(k)}\bigr)\]
Furthermore you can exchange the exponential map, the logarithmic map, and the parallel transport by a retraction, an inverse retraction, and a vector transport.
Finally you can also update the base points $m$ and $n$ during the iterations. This introduces a few additional vector transports. The same holds for the case $Λ(m^{(k)})\neq n^{(k)}$ at some point. All these cases are covered in the algorithm.
Manopt.ChambollePock
— FunctionChambollePock(M, N, f, p, X, m, n, prox_G, prox_G_dual, adjoint_linear_operator; kwargs...)
ChambollePock!(M, N, f, p, X, m, n, prox_G, prox_G_dual, adjoint_linear_operator; kwargs...)
Perform the Riemannian Chambolle—Pock algorithm.
Given a cost
function $\mathcal E:\mathcal M → ℝ$ of the form
\[\mathcal f(p) = F(p) + G( Λ(p) ),\]
where $F:\mathcal M → ℝ$, $G:\mathcal N → ℝ$, and $Λ:\mathcal M → \mathcal N$.
This can be done inplace of $p$.
Input parameters
M::
AbstractManifold
: a Riemannian manifold $\mathcal M$
N::
AbstractManifold
: a Riemannian manifold $\mathcal M$
p
: a point on the manifold $\mathcal M$X
: a tangent vector at the point $p$ on the manifold $\mathcal M$m
: a point on the manifold $\mathcal M$n
: a point on the manifold $\mathcal N$adjoint_linearized_operator
: the adjoint $DΛ^*$ of the linearized operator $DΛ: T_{m}\mathcal M → T_{Λ(m)}\mathcal N)$prox_F, prox_G_Dual
: the proximal maps of $F$ and $G^\ast_n$
note that depending on the AbstractEvaluationType
evaluation
the last three parameters as well as the forward operator Λ
and the linearized_forward_operator
can be given as allocating functions (Manifolds, parameters) -> result
or as mutating functions (Manifold, result, parameters)
-> result` to spare allocations.
By default, this performs the exact Riemannian Chambolle Pock algorithm, see the optional parameter DΛ
for their linearized variant.
For more details on the algorithm, see [BHS+21].
Keyword Arguments
acceleration=0.05
: acceleration parameterdual_stepsize=1/sqrt(8)
: proximal parameter of the primal proxevaluation=
AllocatingEvaluation
()
: specify whether the functions that return an array, for example a point or a tangent vector, work by allocating its result (AllocatingEvaluation
) or whether they modify their input argument to return the result therein (InplaceEvaluation
). Since usually the first argument is the manifold, the modified argument is the second.inverse_retraction_method=
default_inverse_retraction_method
(M, typeof(p))
: an inverse retraction $\operatorname{retr}^{-1}$ to use, see the section on retractions and their inversesinverse_retraction_method_dual=
default_inverse_retraction_method
(N, typeof(n))
: an inverse retraction $\operatorname{retr}^{-1}$ to use, see the section on retractions and their inversesΛ=missing
: the (forward) operator $Λ(⋅)$ (required for the:exact
variant)linearized_forward_operator=missing
: its linearization $DΛ(⋅)[⋅]$ (required for the:linearized
variant)primal_stepsize=1/sqrt(8)
: proximal parameter of the dual proxrelaxation=1.
: the relaxation parameter $γ$relax=:primal
: whether to relax the primal or dualvariant=:exact
ifΛ
is missing, otherwise:linearized
: variant to use. Note that this changes the arguments theforward_operator
is called with.stopping_criterion=
StopAfterIteration`(100)
: a functor indicating that the stopping criterion is fulfilledupdate_primal_base=missing
: function to updatem
(identity by default/missing)update_dual_base=missing
: function to updaten
(identity by default/missing)retraction_method=
default_retraction_method
(M, typeof(p))
: a retraction $\operatorname{retr}$ to use, see the section on retractionsvector_transport_method=
default_vector_transport_method
(M, typeof(p))
: a vector transport $\mathcal T_{⋅←⋅}$ to use, see the section on vector transportsvector_transport_method_dual=
default_vector_transport_method
(N, typeof(n))
: a vector transport $\mathcal T_{⋅←⋅}$ to use, see the section on vector transports
Output
The obtained approximate minimizer $p^*$. To obtain the whole final state of the solver, see get_solver_return
for details, especially the return_state=
keyword.
Manopt.ChambollePock!
— FunctionChambollePock(M, N, f, p, X, m, n, prox_G, prox_G_dual, adjoint_linear_operator; kwargs...)
ChambollePock!(M, N, f, p, X, m, n, prox_G, prox_G_dual, adjoint_linear_operator; kwargs...)
Perform the Riemannian Chambolle—Pock algorithm.
Given a cost
function $\mathcal E:\mathcal M → ℝ$ of the form
\[\mathcal f(p) = F(p) + G( Λ(p) ),\]
where $F:\mathcal M → ℝ$, $G:\mathcal N → ℝ$, and $Λ:\mathcal M → \mathcal N$.
This can be done inplace of $p$.
Input parameters
M::
AbstractManifold
: a Riemannian manifold $\mathcal M$
N::
AbstractManifold
: a Riemannian manifold $\mathcal M$
p
: a point on the manifold $\mathcal M$X
: a tangent vector at the point $p$ on the manifold $\mathcal M$m
: a point on the manifold $\mathcal M$n
: a point on the manifold $\mathcal N$adjoint_linearized_operator
: the adjoint $DΛ^*$ of the linearized operator $DΛ: T_{m}\mathcal M → T_{Λ(m)}\mathcal N)$prox_F, prox_G_Dual
: the proximal maps of $F$ and $G^\ast_n$
note that depending on the AbstractEvaluationType
evaluation
the last three parameters as well as the forward operator Λ
and the linearized_forward_operator
can be given as allocating functions (Manifolds, parameters) -> result
or as mutating functions (Manifold, result, parameters)
-> result` to spare allocations.
By default, this performs the exact Riemannian Chambolle Pock algorithm, see the optional parameter DΛ
for their linearized variant.
For more details on the algorithm, see [BHS+21].
Keyword Arguments
acceleration=0.05
: acceleration parameterdual_stepsize=1/sqrt(8)
: proximal parameter of the primal proxevaluation=
AllocatingEvaluation
()
: specify whether the functions that return an array, for example a point or a tangent vector, work by allocating its result (AllocatingEvaluation
) or whether they modify their input argument to return the result therein (InplaceEvaluation
). Since usually the first argument is the manifold, the modified argument is the second.inverse_retraction_method=
default_inverse_retraction_method
(M, typeof(p))
: an inverse retraction $\operatorname{retr}^{-1}$ to use, see the section on retractions and their inversesinverse_retraction_method_dual=
default_inverse_retraction_method
(N, typeof(n))
: an inverse retraction $\operatorname{retr}^{-1}$ to use, see the section on retractions and their inversesΛ=missing
: the (forward) operator $Λ(⋅)$ (required for the:exact
variant)linearized_forward_operator=missing
: its linearization $DΛ(⋅)[⋅]$ (required for the:linearized
variant)primal_stepsize=1/sqrt(8)
: proximal parameter of the dual proxrelaxation=1.
: the relaxation parameter $γ$relax=:primal
: whether to relax the primal or dualvariant=:exact
ifΛ
is missing, otherwise:linearized
: variant to use. Note that this changes the arguments theforward_operator
is called with.stopping_criterion=
StopAfterIteration`(100)
: a functor indicating that the stopping criterion is fulfilledupdate_primal_base=missing
: function to updatem
(identity by default/missing)update_dual_base=missing
: function to updaten
(identity by default/missing)retraction_method=
default_retraction_method
(M, typeof(p))
: a retraction $\operatorname{retr}$ to use, see the section on retractionsvector_transport_method=
default_vector_transport_method
(M, typeof(p))
: a vector transport $\mathcal T_{⋅←⋅}$ to use, see the section on vector transportsvector_transport_method_dual=
default_vector_transport_method
(N, typeof(n))
: a vector transport $\mathcal T_{⋅←⋅}$ to use, see the section on vector transports
Output
The obtained approximate minimizer $p^*$. To obtain the whole final state of the solver, see get_solver_return
for details, especially the return_state=
keyword.
State
Manopt.ChambollePockState
— TypeChambollePockState <: AbstractPrimalDualSolverState
stores all options and variables within a linearized or exact Chambolle Pock.
Fields
acceleration::R
: acceleration factordual_stepsize::R
: proximal parameter of the dual proxinverse_retraction_method::AbstractInverseRetractionMethod
: an inverse retraction $\operatorname{retr}^{-1}$ to use, see the section on retractions and their inversesinverse_retraction_method_dual::AbstractInverseRetractionMethod
: an inverse retraction $\operatorname{retr}^{-1}$ to use, see the section on retractions and their inversesm::P
: base point on $\mathcal M$n::Q
: base point on $\mathcal N$p::P
: an initial point on $p^{(0)} ∈ \mathcal M$pbar::P
: the relaxed iterate used in the next dual update step (when using:primal
relaxation)primal_stepsize::R
: proximal parameter of the primal proxX::T
: an initial tangent vector $X^{(0)} ∈ T_{p^{(0)}}\mathcal M$Xbar::T
: the relaxed iterate used in the next primal update step (when using:dual
relaxation)relaxation::R
: relaxation in the primal relaxation step (to computepbar
:relax::Symbol: which variable to relax (
:primalor
:dual`:retraction_method::AbstractRetractionMethod
: a retraction $\operatorname{retr}$ to use, see the section on retractionsstop::StoppingCriterion
: a functor indicating that the stopping criterion is fulfilledvariant
: whether to perform an:exact
or:linearized
Chambolle-Pockupdate_primal_base
: function(pr, st, k) -> m
to update the primal baseupdate_dual_base
: function(pr, st, k) -> n
to update the dual basevector_transport_method::AbstractVectorTransportMethodP
: a vector transport $\mathcal T_{⋅←⋅}$ to use, see the section on vector transportsvector_transport_method_dual::AbstractVectorTransportMethodP
: a vector transport $\mathcal T_{⋅←⋅}$ to use, see the section on vector transports
Here, P
is a point type on $\mathcal M$, T
its tangent vector type, Q
a point type on $\mathcal N$, and R<:Real
is a real number type
where for the last two the functions a AbstractManoptProblem
p
, AbstractManoptSolverState
o
and the current iterate i
are the arguments. If you activate these to be different from the default identity, you have to provide p.Λ
for the algorithm to work (which might be missing
in the linearized case).
Constructor
ChambollePockState(M::AbstractManifold, N::AbstractManifold;
kwargs...
) where {P, Q, T, R <: Real}
Keyword arguments
n=
[
rand](@extref Base.rand-Tuple{AbstractManifold})
(N)`p=
rand
(M)
m=
rand
(M)
X=
zero_vector
(M, p)
acceleration=0.0
dual_stepsize=1/sqrt(8)
primal_stepsize=1/sqrt(8)
inverse_retraction_method=
default_inverse_retraction_method
(M, typeof(p))
: an inverse retraction $\operatorname{retr}^{-1}$ to use, see the section on retractions and their inversesinverse_retraction_method_dual=
default_inverse_retraction_method
(N, typeof(n))
: an inverse retraction $\operatorname{retr}^{-1}$ to use, see the section on retractions and their inversesrelaxation=1.0
relax=:primal
: relax the primal variable by defaultretraction_method=
default_retraction_method
(M, typeof(p))
: a retraction $\operatorname{retr}$ to use, see the section on retractionsstopping_criterion=
StopAfterIteration
(300)
: a functor indicating that the stopping criterion is fulfilledvariant=:exact
: run the exact Chambolle Pock by defaultupdate_primal_base=missing
update_dual_base=missing
vector_transport_method=
default_vector_transport_method
(M, typeof(p))
: a vector transport $\mathcal T_{⋅←⋅}$ to use, see the section on vector transportsvector_transport_method_dual=
default_vector_transport_method
(N, typeof(n))
: a vector transport $\mathcal T_{⋅←⋅}$ to use, see the section on vector transports
if Manifolds.jl
is loaded, N
is also a keyword argument and set to TangentBundle(M)
by default.
Useful terms
Manopt.primal_residual
— Functionprimal_residual(p, o, x_old, X_old, n_old)
Compute the primal residual at current iterate $k$ given the necessary values $x_{k-1}, X_{k-1}$, and $n_{k-1}$ from the previous iterate.
\[\Bigl\lVert \frac{1}{σ}\operatorname{retr}^{-1}_{x_{k}}x_{k-1} - V_{x_k\gets m_k}\bigl(DΛ^*(m_k)\bigl[V_{n_k\gets n_{k-1}}X_{k-1} - X_k \bigr] \Bigr\rVert\]
where $V_{⋅\gets⋅}$ is the vector transport used in the ChambollePockState
Manopt.dual_residual
— Functiondual_residual(p, o, x_old, X_old, n_old)
Compute the dual residual at current iterate $k$ given the necessary values $x_{k-1}, X_{k-1}$, and $n_{k-1}$ from the previous iterate. The formula is slightly different depending on the o.variant
used:
For the :linearized
it reads
\[\Bigl\lVert \frac{1}{τ}\bigl( V_{n_{k}\gets n_{k-1}}(X_{k-1}) - X_k \bigr) - DΛ(m_k)\bigl[ V_{m_k\gets x_k}\operatorname{retr}^{-1}_{x_{k}}x_{k-1} \bigr] \Bigr\rVert\]
and for the :exact
variant
\[\Bigl\lVert \frac{1}{τ} V_{n_{k}\gets n_{k-1}}(X_{k-1}) - \operatorname{retr}^{-1}_{n_{k}}\bigl( Λ(\operatorname{retr}_{m_{k}}(V_{m_k\gets x_k}\operatorname{retr}^{-1}_{x_{k}}x_{k-1})) \bigr) \Bigr\rVert\]
where in both cases $V_{⋅\gets⋅}$ is the vector transport used in the ChambollePockState
.
Debug
Manopt.DebugDualBaseIterate
— FunctionDebugDualBaseIterate(io::IO=stdout)
Print the dual base variable by using DebugEntry
, see their constructors for detail. This method is further set display o.n
.
Manopt.DebugDualBaseChange
— FunctionDebugDualChange(; storage=StoreStateAction([:n]), io::IO=stdout)
Print the change of the dual base variable by using DebugEntryChange
, see their constructors for detail, on o.n
.
Manopt.DebugPrimalBaseIterate
— FunctionDebugPrimalBaseIterate()
Print the primal base variable by using DebugEntry
, see their constructors for detail. This method is further set display o.m
.
Manopt.DebugPrimalBaseChange
— FunctionDebugPrimalBaseChange(a::StoreStateAction=StoreStateAction([:m]),io::IO=stdout)
Print the change of the primal base variable by using DebugEntryChange
, see their constructors for detail, on o.n
.
Manopt.DebugDualChange
— TypeDebugDualChange(opts...)
Print the change of the dual variable, similar to DebugChange
, see their constructors for detail, but with a different calculation of the change, since the dual variable lives in (possibly different) tangent spaces.
Manopt.DebugDualIterate
— FunctionDebugDualIterate(e)
Print the dual variable by using DebugEntry
, see their constructors for detail. This method is further set display o.X
.
Manopt.DebugDualResidual
— TypeDebugDualResidual <: DebugAction
A Debug action to print the dual residual. The constructor accepts a printing function and some (shared) storage, which should at least record :Iterate
, :X
and :n
.
Constructor
DebugDualResidual(; kwargs...)
Keyword warguments
io=
stdout`: stream to perform the debug toformat="$prefix%s"
: format to print the dual residual, using theprefix="Dual Residual: "
: short form to just set the prefixstorage
(a newStoreStateAction
) to store values for the debug.
Manopt.DebugPrimalChange
— FunctionDebugPrimalChange(opts...)
Print the change of the primal variable by using DebugChange
, see their constructors for detail.
Manopt.DebugPrimalIterate
— FunctionDebugPrimalIterate(opts...;kwargs...)
Print the change of the primal variable by using DebugIterate
, see their constructors for detail.
Manopt.DebugPrimalResidual
— TypeDebugPrimalResidual <: DebugAction
A Debug action to print the primal residual. The constructor accepts a printing function and some (shared) storage, which should at least record :Iterate
, :X
and :n
.
Constructor
DebugPrimalResidual(; kwargs...)
Keyword warguments
io=
stdout`: stream to perform the debug toformat="$prefix%s"
: format to print the dual residual, using theprefix="Primal Residual: "
: short form to just set the prefixstorage
(a newStoreStateAction
) to store values for the debug.
Manopt.DebugPrimalDualResidual
— TypeDebugPrimalDualResidual <: DebugAction
A Debug action to print the primal dual residual. The constructor accepts a printing function and some (shared) storage, which should at least record :Iterate
, :X
and :n
.
Constructor
DebugPrimalDualResidual()
with the keywords
Keyword warguments
io=
stdout`: stream to perform the debug toformat="$prefix%s"
: format to print the dual residual, using theprefix="PD Residual: "
: short form to just set the prefixstorage
(a newStoreStateAction
) to store values for the debug.
Record
Manopt.RecordDualBaseIterate
— FunctionRecordDualBaseIterate(n)
Create an RecordAction
that records the dual base point, an RecordEntry
of o.n
.
Manopt.RecordDualBaseChange
— FunctionRecordDualBaseChange(e)
Create an RecordAction
that records the dual base point change, an RecordEntryChange
of o.n
with distance to the last value to store a value.
Manopt.RecordDualChange
— FunctionRecordDualChange()
Create the action either with a given (shared) Storage, which can be set to the values
Tuple, if that is provided).
Manopt.RecordDualIterate
— FunctionRecordDualIterate(X)
Create an RecordAction
that records the dual base point, an RecordEntry
of o.X
.
Manopt.RecordPrimalBaseIterate
— FunctionRecordPrimalBaseIterate(x)
Create an RecordAction
that records the primal base point, an RecordEntry
of o.m
.
Manopt.RecordPrimalBaseChange
— FunctionRecordPrimalBaseChange()
Create an RecordAction
that records the primal base point change, an RecordEntryChange
of o.m
with distance to the last value to store a value.
Manopt.RecordPrimalChange
— FunctionRecordPrimalChange(a)
Create an RecordAction
that records the primal value change, RecordChange
, to record the change of o.x
.
Manopt.RecordPrimalIterate
— FunctionRecordDualBaseIterate(x)
Create an RecordAction
that records the dual base point, an RecordIterate
of o.x
.
Internals
Manopt.update_prox_parameters!
— Functionupdate_prox_parameters!(o)
update the prox parameters as described in Algorithm 2 of [CP11],
- $θ_{n} = \frac{1}{\sqrt{1+2γτ_n}}$
- $τ_{n+1} = θ_nτ_n$
- $σ_{n+1} = \frac{σ_n}{θ_n}$
Technical details
The ChambollePock
solver requires the following functions of a manifold to be available for both the manifold $\mathcal M$and $\mathcal N$
- A
retract!
(M, q, p, X)
; it is recommended to set thedefault_retraction_method
to a favourite retraction. If this default is set, aretraction_method=
orretraction_method_dual=
(for $\mathcal N$) does not have to be specified. - An
inverse_retract!
(M, X, p, q)
; it is recommended to set thedefault_inverse_retraction_method
to a favourite retraction. If this default is set, ainverse_retraction_method=
orinverse_retraction_method_dual=
(for $\mathcal N$) does not have to be specified. - A
vector_transport_to!
M, Y, p, X, q)
; it is recommended to set thedefault_vector_transport_method
to a favourite retraction. If this default is set, avector_transport_method=
orvector_transport_method_dual=
(for $\mathcal N$) does not have to be specified. - A `copyto!
(M, q, p)
andcopy
(M,p)
for points.
Literature
- [BHS+21]
- R. Bergmann, R. Herzog, M. Silva Louzeiro, D. Tenbrinck and J. Vidal-Núñez. Fenchel duality theory and a primal-dual algorithm on Riemannian manifolds. Foundations of Computational Mathematics 21, 1465–1504 (2021), arXiv:1908.02022.
- [CP11]
- A. Chambolle and T. Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision 40, 120–145 (2011).