diff --git a/cpp/sopt/CMakeLists.txt b/cpp/sopt/CMakeLists.txt index 4b01b5bb..69bdf71d 100644 --- a/cpp/sopt/CMakeLists.txt +++ b/cpp/sopt/CMakeLists.txt @@ -8,6 +8,7 @@ set(headers maths.h proximal.h relative_variation.h sdmm.h wavelets.h conjugate_gradient.h l1_proximal.h logging.enabled.h padmm.h proximal_expression.h reweighted.h types.h wrapper.h exception.h linear_transform.h logging.h positive_quadrant.h + l2_forward_backward.h l2_primal_dual.h real_type.h sampling.h power_method.h objective_functions.h ${PROJECT_BINARY_DIR}/include/sopt/config.h ) set(wavelet_headers diff --git a/cpp/sopt/l2_forward_backward.h b/cpp/sopt/l2_forward_backward.h new file mode 100644 index 00000000..3786b315 --- /dev/null +++ b/cpp/sopt/l2_forward_backward.h @@ -0,0 +1,341 @@ +#ifndef SOPT_L2_FORWARD_BACKWARD_H +#define SOPT_L2_FORWARD_BACKWARD_H + +#include "sopt/config.h" +#include +#include +#include +#include "sopt/exception.h" +#include "sopt/forward_backward.h" +#include "sopt/l1_proximal.h" +#include "sopt/linear_transform.h" +#include "sopt/logging.h" +#include "sopt/proximal.h" +#include "sopt/relative_variation.h" +#include "sopt/types.h" + +#ifdef SOPT_MPI +#include "sopt/mpi/communicator.h" +#include "sopt/mpi/utilities.h" +#endif + +namespace sopt { +namespace algorithm { +template +class ImagingForwardBackward { + //! Underlying algorithm + typedef ForwardBackward FB; + + public: + typedef typename FB::value_type value_type; + typedef typename FB::Scalar Scalar; + typedef typename FB::Real Real; + typedef typename FB::t_Vector t_Vector; + typedef typename FB::t_LinearTransform t_LinearTransform; + template + using t_Proximal = std::function; + typedef typename FB::t_Gradient t_Gradient; + typedef typename FB::t_IsConverged t_IsConverged; + + //! Values indicating how the algorithm ran + struct Diagnostic : public FB::Diagnostic { + Diagnostic(t_uint niters = 0u, bool good = false) : FB::Diagnostic(niters, good) {} + Diagnostic(t_uint niters, bool good, t_Vector &&residual) + : PD::Diagnostic(niters, good, std::move(residual)) {} + }; + //! Holds result vector as well + struct DiagnosticAndResult : public Diagnostic { + //! Output x + t_Vector x; + }; + + //! Setups imaging wrapper for ForwardBackward + //! \param[in] f_proximal: proximal operator of the \f$f\f$ function. + //! \param[in] g_proximal: proximal operator of the \f$g\f$ function + template + ImagingForwardBackward(Eigen::MatrixBase const &target) + : l2_proximal_([](t_Vector &output, const t_real &gamma, const t_Vector &x) -> void { + proximal::l2_norm(output, gamma, x); + }), + l2_proximal_weighted_( + [](t_Vector &output, const Vector &gamma, const t_Vector &x) -> void { + proximal::l2_norm(output, gamma, x); + }), + l2_proximal_weights_(Vector::Ones(1)), + l2_gradient_([](t_Vector &output, const t_Vector &x) -> void { + output = x; + }), // gradient of 1/2 * x^2 = x; + tight_frame_(false), + residual_tolerance_(0.), + relative_variation_(1e-4), + residual_convergence_(nullptr), + objective_convergence_(nullptr), + itermax_(std::numeric_limits::max()), + gamma_(1e-8), + beta_(1), + sigma_(1), + nu_(1), + is_converged_(), + Phi_(linear_transform_identity()), + target_(target) {} + virtual ~ImagingForwardBackward() {} + +// Macro helps define properties that can be initialized as in +// auto padmm = ImagingForwardBackward().prop0(value).prop1(value); +#define SOPT_MACRO(NAME, TYPE) \ + TYPE const &NAME() const { return NAME##_; } \ + ImagingForwardBackward &NAME(TYPE const &NAME) { \ + NAME##_ = NAME; \ + return *this; \ + } \ + \ + protected: \ + TYPE NAME##_; \ + \ + public: + + //! l2 proximal for regularizaiton + SOPT_MACRO(l2_proximal, t_Proximal); + //! l2 proximal for regularizaiton with weights + SOPT_MACRO(l2_proximal, t_Proximal>); + //! l2 proximal weights + SOPT_MACRO(l2_proximal_weights, Vector); + //! Gradient of the l2 norm + SOPT_MACRO(l2_gradient, t_Gradient); + //! Whether Ψ is a tight-frame or not + SOPT_MACRO(tight_frame, bool); + //! \brief Convergence of the relative variation of the objective functions + //! \details If negative, this convergence criteria is disabled. + SOPT_MACRO(residual_tolerance, Real); + //! \brief Convergence of the relative variation of the objective functions + //! \details If negative, this convergence criteria is disabled. + SOPT_MACRO(relative_variation, Real); + //! \brief Convergence of the residuals + //! \details If negative, this convergence criteria is disabled. + SOPT_MACRO(residual_convergence, t_IsConverged); + //! \brief Convergence of the residuals + //! \details If negative, this convergence criteria is disabled. + SOPT_MACRO(objective_convergence, t_IsConverged); + //! Maximum number of iterations + SOPT_MACRO(itermax, t_uint); + //! γ parameter + SOPT_MACRO(gamma, Real); + //! γ parameter + SOPT_MACRO(beta, Real); + //! γ parameter + SOPT_MACRO(sigma, Real); + //! ν parameter + SOPT_MACRO(nu, Real); + //! A function verifying convergence + SOPT_MACRO(is_converged, t_IsConverged); + //! Measurement operator + SOPT_MACRO(Phi, t_LinearTransform); +#ifdef SOPT_MPI + //! Communicator for summing objective_function + SOPT_MACRO(obj_comm, mpi::Communicator); +#endif + +#undef SOPT_MACRO + //! Vector of target measurements + t_Vector const &target() const { return target_; } + //! Minimun of objective_function + Real objmin() const { return objmin_; }; + //! Sets the vector of target measurements + template + ImagingForwardBackward &target(Eigen::MatrixBase const &target) { + target_ = target; + return *this; + } + + //! \brief Calls Forward Backward + //! \param[out] out: Output vector x + Diagnostic operator()(t_Vector &out) const { + return operator()(out, ForwardBackward::initial_guess(target(), Phi(), nu())); + } + //! \brief Calls Forward Backward + //! \param[out] out: Output vector x + //! \param[in] guess: initial guess + Diagnostic operator()(t_Vector &out, std::tuple const &guess) const { + return operator()(out, std::get<0>(guess), std::get<1>(guess)); + } + //! \brief Calls Forward Backward + //! \param[out] out: Output vector x + //! \param[in] guess: initial guess + Diagnostic operator()(t_Vector &out, + std::tuple const &guess) const { + return operator()(out, std::get<0>(guess), std::get<1>(guess)); + } + //! \brief Calls Forward Backward + //! \param[in] guess: initial guess + DiagnosticAndResult operator()(std::tuple const &guess) const { + return operator()(std::tie(std::get<0>(guess), std::get<1>(guess))); + } + //! \brief Calls Forward Backward + //! \param[in] guess: initial guess + DiagnosticAndResult operator()( + std::tuple const &guess) const { + DiagnosticAndResult result; + static_cast(result) = operator()(result.x, guess); + return result; + } + //! \brief Calls Forward Backward + //! \param[in] guess: initial guess + DiagnosticAndResult operator()() const { + DiagnosticAndResult result; + static_cast(result) = operator()( + result.x, ForwardBackward::initial_guess(target(), Phi(), nu())); + return result; + } + //! Makes it simple to chain different calls to FB + DiagnosticAndResult operator()(DiagnosticAndResult const &warmstart) const { + DiagnosticAndResult result = warmstart; + static_cast(result) = operator()(result.x, warmstart.x, warmstart.residual); + return result; + } + + //! Set Φ and Φ^† using arguments that sopt::linear_transform understands + template + typename std::enable_if= 1, ImagingForwardBackward &>::type Phi( + ARGS &&... args) { + Phi_ = linear_transform(std::forward(args)...); + return *this; + } + + //! \brief L1 proximal used during calculation + //! \details Non-const version to setup the object. + t_Proximal &l1_proximal() { return l1_proximal_; } + //! \brief Proximal of the L2 ball + //! \details Non-const version to setup the object. + t_Gradient &l2_graident() { return l2_gradient_; } + + //! Helper function to set-up default residual convergence function + ImagingForwardBackward &residual_convergence(Real const &tolerance) { + return residual_convergence(nullptr).residual_tolerance(tolerance); + } + //! Helper function to set-up default residual convergence function + ImagingForwardBackward &objective_convergence(Real const &tolerance) { + return objective_convergence(nullptr).relative_variation(tolerance); + } + //! Convergence function that takes only the output as argument + ImagingForwardBackward &is_converged(std::function const &func) { + return is_converged([func](t_Vector const &x, t_Vector const &) { return func(x); }); + } + + protected: + //! Vector of measurements + t_Vector target_; + //! Mininum of objective function + mutable Real objmin_; + + //! \brief Calls Forward Backward + //! \param[out] out: Output vector x + //! \param[in] guess: initial guess + //! \param[in] residuals: initial residuals + Diagnostic operator()(t_Vector &out, t_Vector const &guess, t_Vector const &res) const; + + //! Helper function to simplify checking convergence + bool residual_convergence(t_Vector const &x, t_Vector const &residual) const; + + //! Helper function to simplify checking convergence + bool objective_convergence(ScalarRelativeVariation &scalvar, t_Vector const &x, + t_Vector const &residual) const; +#ifdef SOPT_MPI + //! Helper function to simplify checking convergence + bool objective_convergence(mpi::Communicator const &obj_comm, + ScalarRelativeVariation &scalvar, t_Vector const &x, + t_Vector const &residual) const; +#endif + + //! Helper function to simplify checking convergence + bool is_converged(ScalarRelativeVariation &scalvar, t_Vector const &x, + t_Vector const &residual) const; +}; + +template +typename ImagingForwardBackward::Diagnostic ImagingForwardBackward::operator()( + t_Vector &out, t_Vector const &guess, t_Vector const &res) const { + SOPT_HIGH_LOG("Performing Forward Backward with L2 and L2 norms"); + // The f proximal is an L2 proximal + auto const g_proximal = [this](t_Vector &out, Real gamma, t_Vector const &x) { + if (this->l2_proximal_weights().size() > 1) + this->l2_proximal_weighted()(out, this->l1_proximal_weights() * gamma, x); + else + this->l2_proximal()(out, this->l1_proximal_weights()(0) * gamma, x); + }; + const Real sigma_factor = sigma() * sigma(); + auto const f_gradient = [this, sigma_factor](t_Vector &out, t_Vector const &x) { + this->l2_gradient()(out, x / sigma_factor); + }; + ScalarRelativeVariation scalvar(relative_variation(), relative_variation(), + "Objective function"); + auto const convergence = [this, &scalvar](t_Vector const &x, t_Vector const &residual) mutable { + const bool result = this->is_converged(scalvar, x, residual); + this->objmin_ = std::real(scalvar.previous()); + return result; + }; + auto const fb = ForwardBackward(f_gradient, g_proximal, target()) + .itermax(itermax()) + .beta(beta()) + .gamma(gamma()) + .nu(nu()) + .Phi(Phi()) + .is_converged(convergence); + static_cast::Diagnostic &>(result) = + fb(out, std::tie(guess, res)); + return result; +} + +template +bool ImagingForwardBackward::residual_convergence(t_Vector const &x, + t_Vector const &residual) const { + if (static_cast(residual_convergence())) return residual_convergence()(x, residual); + if (residual_tolerance() <= 0e0) return true; + auto const residual_norm = sopt::l2_norm(residual); + SOPT_LOW_LOG(" - [FB] Residuals: {} +bool ImagingForwardBackward::objective_convergence(ScalarRelativeVariation &scalvar, + t_Vector const &x, + t_Vector const &residual) const { + if (static_cast(objective_convergence())) return objective_convergence()(x, residual); + if (scalvar.relative_tolerance() <= 0e0) return true; + auto const current = ((gamma() > 0) ? sopt::l2_norm(x, l2_proximal_weights()) * gamma() : 0) + + std::pow(sopt::l2_norm(residual), 2) / (2 * sigma() * sigma()); + return scalvar(current); +}; + +#ifdef SOPT_MPI +template +bool ImagingForwardBackward::objective_convergence(mpi::Communicator const &obj_comm, + ScalarRelativeVariation &scalvar, + t_Vector const &x, + t_Vector const &residual) const { + if (static_cast(objective_convergence())) return objective_convergence()(x, residual); + if (scalvar.relative_tolerance() <= 0e0) return true; + auto const current = obj_comm.all_sum_all( + ((gamma() > 0) ? sopt::l2_norm(x, l2_proximal_weights()) * gamma() : 0) + + std::pow(sopt::l2_norm(residual), 2) / (2 * sigma() * sigma())); + return scalvar(current); +}; +#endif + +template +bool ImagingForwardBackward::is_converged(ScalarRelativeVariation &scalvar, + t_Vector const &x, + t_Vector const &residual) const { + auto const user = static_cast(is_converged()) == false or is_converged()(x, residual); + auto const res = residual_convergence(x, residual); +#ifdef SOPT_MPI + auto const obj = objective_convergence(obj_comm(), scalvar, x, residual); +#else + auto const obj = objective_convergence(scalvar, x, residual); +#endif + // beware of short-circuiting! + // better evaluate each convergence function everytime, especially with mpi + return user and res and obj; +} +} // namespace algorithm +} // namespace sopt +#endif diff --git a/cpp/sopt/l2_primal_dual.h b/cpp/sopt/l2_primal_dual.h new file mode 100644 index 00000000..39670f52 --- /dev/null +++ b/cpp/sopt/l2_primal_dual.h @@ -0,0 +1,400 @@ +#ifndef SOPT_L2_PRIMAL_DUAL_H +#define SOPT_L2_PRIMAL_DUAL_H + +#include "sopt/config.h" +#include +#include +#include +#include "sopt/exception.h" +#include "sopt/linear_transform.h" +#include "sopt/logging.h" +#include "sopt/primal_dual.h" +#include "sopt/proximal.h" +#include "sopt/relative_variation.h" +#include "sopt/types.h" + +namespace sopt { +namespace algorithm { +template +class ImagingPrimalDual { + //! Underlying algorithm + typedef PrimalDual PD; + + public: + typedef typename PD::value_type value_type; + typedef typename PD::Scalar Scalar; + typedef typename PD::Real Real; + typedef typename PD::t_Vector t_Vector; + typedef typename PD::t_LinearTransform t_LinearTransform; + template + using t_Proximal = std::function; + typedef typename PD::t_IsConverged t_IsConverged; + typedef typename PD::t_Constraint t_Constraint; + typedef typename PD::t_Random_Updater t_Random_Updater; + + //! Values indicating how the algorithm ran + struct Diagnostic : public PD::Diagnostic { + Diagnostic(t_uint niters = 0u, bool good = false) : PD::Diagnostic(niters, good) {} + Diagnostic(t_uint niters, bool good, t_Vector &&residual) + : PD::Diagnostic(niters, good, std::move(residual)) {} + }; + //! Holds result vector as well + struct DiagnosticAndResult : public Diagnostic { + //! Output x + t_Vector x; + }; + + //! Setups imaging wrapper for PD + //! \param[in] f_proximal: proximal operator of the \f$f\f$ function. + //! \param[in] g_proximal: proximal operator of the \f$g\f$ function + template + ImagingPrimalDual(Eigen::MatrixBase const &target) + : l2_proximal_([](t_Vector &out, const Real &gamma, const t_Vector &x) { + proximal::l2_norm(out, gamma, x); + }), + l2_proximal_weighted_([](t_Vector &out, const Vector &gamma, const t_Vector &x) { + proximal::l2_norm(out, gamma, x); + }), + l2_proximal_weights_(Vector::Ones(1)), + l2ball_proximal_(1e0), + residual_tolerance_(1e-4), + relative_variation_(1e-4), + residual_convergence_(nullptr), + objective_convergence_(nullptr), + itermax_(std::numeric_limits::max()), + sigma_(1), + tau_(0.5), + gamma_(0.5), + update_scale_(1), + precondition_stepsize_(0.5), + precondition_weights_(t_Vector::Ones(target.size())), + precondition_iters_(0), + xi_(1), + rho_(1), + nu_(1), + is_converged_(), + Phi_(linear_transform_identity()), + Psi_(linear_transform_identity()), + random_measurement_updater_([]() { return true; }), + random_wavelet_updater_([]() { return true; }), + target_(target) {} + virtual ~ImagingPrimalDual() {} + +// Macro helps define properties that can be initialized as in +// auto padmm = ImagingPrimalDual().prop0(value).prop1(value); +#define SOPT_MACRO(NAME, TYPE) \ + TYPE const &NAME() const { return NAME##_; } \ + ImagingPrimalDual &NAME(TYPE const &NAME) { \ + NAME##_ = NAME; \ + return *this; \ + } \ + \ + protected: \ + TYPE NAME##_; \ + \ + public: + //! The l2 prox functioning as f + SOPT_MACRO(l2_proximal, t_Proximal); + //! The l2 prox with weights functioning as f + SOPT_MACRO(l2_proximal_weighted, t_Proximal>); + //! The l2 prox weights functioning + SOPT_MACRO(l2_proximal_weights, Vector); + //! The weighted L2 proximal functioning as g + SOPT_MACRO(l2ball_proximal, proximal::WeightedL2Ball); + //! \brief Convergence of the relative variation of the objective functions + //! \details If negative, this convergence criteria is disabled. + SOPT_MACRO(residual_tolerance, Real); + //! \brief Convergence of the relative variation of the objective functions + //! \details If negative, this convergence criteria is disabled. + SOPT_MACRO(relative_variation, Real); + //! \brief Convergence of the residuals + //! \details If negative, this convergence criteria is disabled. + SOPT_MACRO(residual_convergence, t_IsConverged); + //! \brief Convergence of the residuals + //! \details If negative, this convergence criteria is disabled. + SOPT_MACRO(objective_convergence, t_IsConverged); + //! Maximum number of iterations + SOPT_MACRO(itermax, t_uint); + //! gamma parameter + SOPT_MACRO(gamma, Real); + //! update parameter + SOPT_MACRO(update_scale, Real); + //! Apply positivity constraint + SOPT_MACRO(positivity_constraint, bool); + //! Apply real constraint + SOPT_MACRO(real_constraint, bool); + //! sigma parameter + SOPT_MACRO(sigma, Real); + //! tau parameter + SOPT_MACRO(tau, Real); + //! xi parameter + SOPT_MACRO(xi, Real); + //! rho parameter + SOPT_MACRO(rho, Real); + //! ν parameter + SOPT_MACRO(nu, Real); + //! precondtion step size parameter + SOPT_MACRO(precondition_stepsize, Real); + //! precondition weights parameter + SOPT_MACRO(precondition_weights, t_Vector); + //! precondition iterations parameter + SOPT_MACRO(precondition_iters, t_uint); + //! A function verifying convergence + SOPT_MACRO(is_converged, t_IsConverged); + //! Measurement operator + SOPT_MACRO(Phi, t_LinearTransform); + //! Wavelet operator + SOPT_MACRO(Psi, t_LinearTransform); + //! lambda that determines if to update measurements + SOPT_MACRO(random_measurement_updater, t_Random_Updater); + //! lambda that determines if to update wavelets + SOPT_MACRO(random_wavelet_updater, t_Random_Updater); +#ifdef SOPT_MPI + //! + SOPT_MACRO(v_all_sum_all_comm, mpi::Communicator); + //! + SOPT_MACRO(u_all_sum_all_comm, mpi::Communicator); +#endif + +#undef SOPT_MACRO + //! Vector of target measurements + t_Vector const &target() const { return target_; } + //! Sets the vector of target measurements + template + ImagingPrimalDual &target(Eigen::MatrixBase const &target) { + target_ = target; + return *this; + } + + //! \brief Calls Primal Dual + //! \param[out] out: Output vector x + Diagnostic operator()(t_Vector &out) const { + return operator()(out, PD::initial_guess(target(), Phi(), nu())); + } + //! \brief Calls Primal Dual + //! \param[out] out: Output vector x + //! \param[in] guess: initial guess + Diagnostic operator()(t_Vector &out, std::tuple const &guess) const { + return operator()(out, std::get<0>(guess), std::get<1>(guess)); + } + //! \brief Calls Primal Dual + //! \param[out] out: Output vector x + //! \param[in] guess: initial guess + Diagnostic operator()(t_Vector &out, + std::tuple const &guess) const { + return operator()(out, std::get<0>(guess), std::get<1>(guess)); + } + //! \brief Calls Primal Dual + //! \param[in] guess: initial guess + DiagnosticAndResult operator()(std::tuple const &guess) const { + return operator()(std::tie(std::get<0>(guess), std::get<1>(guess))); + } + //! \brief Calls Primal Dual + //! \param[in] guess: initial guess + DiagnosticAndResult operator()( + std::tuple const &guess) const { + DiagnosticAndResult result; + static_cast(result) = operator()(result.x, guess); + return result; + } + //! \brief Calls Primal Dual + //! \param[in] guess: initial guess + DiagnosticAndResult operator()() const { + DiagnosticAndResult result; + static_cast(result) = operator()(result.x, + PD::initial_guess(target(), Phi(), nu())); + return result; + } + //! Makes it simple to chain different calls to PD + DiagnosticAndResult operator()(DiagnosticAndResult const &warmstart) const { + DiagnosticAndResult result = warmstart; + static_cast(result) = operator()(result.x, warmstart.x, warmstart.residual); + return result; + } + + //! Set Φ and Φ^† using arguments that sopt::linear_transform understands + template + typename std::enable_if= 1, ImagingPrimalDual &>::type Phi(ARGS &&... args) { + Phi_ = linear_transform(std::forward(args)...); + return *this; + } + + //! \brief Proximal of the L2 ball + //! \details Non-const version to setup the object. + proximal::WeightedL2Ball &l2ball_proximal() { return l2ball_proximal_; } + + //! Analysis operator Ψ + template + typename std::enable_if= 1, ImagingPrimalDual &>::type Psi(ARGS &&... args) { + Psi_ = linear_transform(std::forward(args)...); + return *this; + } + +// Forwards get/setters to L2 and L2Ball proximals +// In practice, we end up with a bunch of functions that make it simpler to set or get values +// associated with the two proximal operators. +// E.g.: `paddm.l2_proximal_itermax(100).l2ball_epsilon(1e-2).l2_proximal_tolerance(1e-4)`. +// ~~~ +#define SOPT_MACRO(VAR, NAME, PROXIMAL) \ + /** \brief Forwards to l2ball_proximal **/ \ + decltype(std::declval const>().VAR()) NAME##_proximal_##VAR() const { \ + return NAME##_proximal().VAR(); \ + } \ + /** \brief Forwards to l2ball_proximal **/ \ + ImagingPrimalDual &NAME##_proximal_##VAR( \ + decltype(std::declval const>().VAR()) VAR) { \ + NAME##_proximal().VAR(VAR); \ + return *this; \ + } + SOPT_MACRO(epsilon, l2ball, WeightedL2Ball); + SOPT_MACRO(weights, l2ball, WeightedL2Ball); +#ifdef SOPT_MPI + SOPT_MACRO(communicator, l2ball, WeightedL2Ball); +#endif +#undef SOPT_MACRO + + //! Helper function to set-up default residual convergence function + ImagingPrimalDual &residual_convergence(Real const &tolerance) { + return residual_convergence(nullptr).residual_tolerance(tolerance); + } + //! Helper function to set-up default residual convergence function + ImagingPrimalDual &objective_convergence(Real const &tolerance) { + return objective_convergence(nullptr).relative_variation(tolerance); + } + //! Convergence function that takes only the output as argument + ImagingPrimalDual &is_converged(std::function const &func) { + return is_converged([func](t_Vector const &x, t_Vector const &) { return func(x); }); + } + + protected: + //! Vector of measurements + t_Vector target_; + + //! \brief Calls Primal Dual + //! \param[out] out: Output vector x + //! \param[in] guess: initial guess + //! \param[in] residuals: initial residuals + Diagnostic operator()(t_Vector &out, t_Vector const &guess, t_Vector const &res) const; + + //! Helper function to simplify checking convergence + bool residual_convergence(t_Vector const &x, t_Vector const &residual) const; + + //! Helper function to simplify checking convergence + bool objective_convergence(ScalarRelativeVariation &scalvar, t_Vector const &x, + t_Vector const &residual) const; + + //! Helper function to simplify checking convergence + bool is_converged(ScalarRelativeVariation &scalvar, t_Vector const &x, + t_Vector const &residual) const; + //! check that l1 and weighted l1 proximal operators are the same function (except for weights) + bool check_l2_weight_proximal(const t_Proximal &no_weights, + const t_Proximal> &with_weights) const { + Vector output; + Vector outputw; + + const Vector x = Vector::Ones(this->l2_proximal_weights().size()); + no_weights(output, 1, x); + with_weights(outputw, Vector::Ones(1), x); + return output.isApprox(outputw); + }; +}; + +template +typename ImagingPrimalDual::Diagnostic ImagingPrimalDual::operator()( + t_Vector &out, t_Vector const &guess, t_Vector const &res) const { + SOPT_HIGH_LOG("Performing Primal Dual with L1 and L2 operators"); + // The f proximal is an L1 proximal that stores some diagnostic result + if (not check_l2_weight_proximal(l2_proximal(), l2_proximal_weighted())) + SOPT_THROW( + "l1 proximal and weighted l1 proximal appear to be different functions. Please make sure " + "both are the same function."); + auto const f_proximal = [this](t_Vector &out, Real gamma, t_Vector const &x) { + if (this->l2_proximal_weights().size() > 1) + this->l2_proximal_weighted()(out, this->l2_proximal_weights() * gamma, x); + else + this->l2_proximal()(out, this->l2_proximal_weights()(0) * gamma, x); + }; + auto const g_proximal = [this](t_Vector &out, Real gamma, t_Vector const &x) { + this->l2ball_proximal()(out, gamma, x); + // applying preconditioning + for (t_int i = 0; i < this->precondition_iters(); i++) + this->l2ball_proximal()( + out, gamma, + out - this->precondition_stepsize() * + (out.array() * this->precondition_weights().array() - x.array()).matrix()); + + if (this->precondition_iters() > 0) out = out.array() * this->precondition_weights().array(); + }; + ScalarRelativeVariation scalvar(relative_variation(), relative_variation(), + "Objective function"); + auto const convergence = [this, scalvar](t_Vector const &x, t_Vector const &residual) mutable { + return this->is_converged(scalvar, x, residual); + }; + const bool positive = positivity_constraint(); + const bool real = real_constraint(); + t_Constraint constraint = [real, positive](t_Vector &out, const t_Vector &x) { + if (real) out.real() = x.real(); + if (positive) out = sopt::positive_quadrant(x); + if (not real and not positive) out = x; + }; + auto const pd = PD(f_proximal, g_proximal, target()) + .itermax(itermax()) + .constraint(constraint) + .sigma(sigma()) + .tau(tau()) + .gamma(gamma()) + .update_scale(update_scale()) + .xi(xi()) + .rho(rho()) + .nu(nu()) + .gamma(gamma()) + .Phi(Phi()) + .Psi(Psi()) + .random_measurement_updater(random_measurement_updater()) + .random_wavelet_updater(random_wavelet_updater()) +#ifdef SOPT_MPI + .v_all_sum_all_comm(v_all_sum_all_comm()) + .u_all_sum_all_comm(u_all_sum_all_comm()) +#endif + .is_converged(convergence); + Diagnostic result; + static_cast(result) = pd(out, std::tie(guess, res)); + return result; +} + +template +bool ImagingPrimalDual::residual_convergence(t_Vector const &x, + t_Vector const &residual) const { + if (static_cast(residual_convergence())) return residual_convergence()(x, residual); + if (residual_tolerance() <= 0e0) return true; + auto const residual_norm = sopt::l2_norm(residual, l2ball_proximal_weights()); + SOPT_LOW_LOG(" - [Primal Dual] Residuals: {} +bool ImagingPrimalDual::objective_convergence(ScalarRelativeVariation &scalvar, + t_Vector const &x, + t_Vector const &residual) const { + if (static_cast(objective_convergence())) return objective_convergence()(x, residual); + if (scalvar.relative_tolerance() <= 0e0) return true; + auto const current = + (l2_proximal_weights().size() > 1) + ? sopt::l2_norm(l2_proximal_weights().array() * (Psi().adjoint() * x).array()) + : sopt::l2_norm(l2_proximal_weights()(0) * (Psi().adjoint() * x)); + return scalvar(current); +}; + +template +bool ImagingPrimalDual::is_converged(ScalarRelativeVariation &scalvar, + t_Vector const &x, t_Vector const &residual) const { + auto const user = static_cast(is_converged()) == false or is_converged()(x, residual); + auto const res = residual_convergence(x, residual); + auto const obj = objective_convergence(scalvar, x, residual); + // beware of short-circuiting! + // better evaluate each convergence function everytime, especially with mpi + return user and res and obj; +} +} // namespace algorithm +} // namespace sopt +#endif diff --git a/cpp/sopt/proximal.h b/cpp/sopt/proximal.h index 9aa9284c..e7d28e78 100644 --- a/cpp/sopt/proximal.h +++ b/cpp/sopt/proximal.h @@ -87,6 +87,12 @@ void l2_norm(Eigen::DenseBase &out, typename real_type: out = x.derived() * 1. / (1. + gamma); } +template +void l2_norm(Eigen::DenseBase &out, Eigen::DenseBase const &gamma, + Eigen::DenseBase const &x) { + out = x.derived().array() * 1. / (1. + gamma).array(); +} + //! Proximal of a function that is always zero, the identity template void id(Eigen::DenseBase &out, typename real_type::type gamma,