proximal.h

Go to the documentation of this file.

#ifndef SOPT_PROXIMAL_H
#define SOPT_PROXIMAL_H
 
#include "sopt/config.h"
#include <type_traits>
#include <Eigen/Core>
#include "sopt/maths.h"
#include "sopt/proximal_expression.h"
#ifdef SOPT_MPI
#include "sopt/mpi/communicator.h"
#include "sopt/mpi/utilities.h"
#endif
 
namespace sopt::proximal {
 
class EuclidianNorm {
 public:
#ifdef SOPT_MPI
  EuclidianNorm(mpi::Communicator const &comm = mpi::Communicator()) : comm_(comm){};
  mpi::Communicator communicator() const { return comm_; }
  EuclidianNorm &communicator(mpi::Communicator const &comm) {
    comm_ = comm;
    return *this;
  }
#endif
  template <typename T0>
  void operator()(Vector<typename T0::Scalar> &out,
                  typename real_type<typename T0::Scalar>::type const &t,
                  Eigen::MatrixBase<T0> const &x) const {
    using Scalar = typename T0::Scalar;
#ifdef SOPT_MPI
    auto const norm = mpi::l2_norm(x, comm_);
#else
    auto const norm = x.norm();
#endif
    if (norm > t)
      out = (Scalar(1) - t / norm) * x;
    else
      out.fill(0);
  }
  template <typename T0>
  ProximalExpression<EuclidianNorm, T0> operator()(typename T0::Scalar const &t,
                                                   Eigen::MatrixBase<T0> const &x) const {
    return {*this, t, x};
  }
#ifdef SOPT_MPI
 private:
  mpi::Communicator comm_;
#endif
};
 
template <typename T0>
auto euclidian_norm(typename real_type<typename T0::Scalar>::type const &t,
                    Eigen::MatrixBase<T0> const &x) -> decltype(EuclidianNorm(), t, x) {
  return EuclidianNorm()(t, x);
}
 
template <typename T0, typename T1>
void l1_norm(Eigen::DenseBase<T0> &out, typename real_type<typename T0::Scalar>::type gamma,
             Eigen::DenseBase<T1> const &x) {
  out = sopt::soft_threshhold<T0>(x, gamma);
}
template <typename T0, typename T1, typename T2>
void l1_norm(Eigen::DenseBase<T0> &out, Eigen::DenseBase<T2> const &gamma,
             Eigen::DenseBase<T1> const &x) {
  out = sopt::soft_threshhold<T0, T2>(x, gamma);
}
 
template <typename S>
void l1_norm(Vector<S> &out, typename real_type<S>::type gamma, Vector<S> const &x) {
  l1_norm<Vector<S>, Vector<S>>(out, gamma, x);
}
 
template <typename T0, typename T1>
void l2_norm(Eigen::DenseBase<T0> &out, typename real_type<typename T0::Scalar>::type gamma,
             Eigen::DenseBase<T1> const &x) {
  out = x.derived() * 1. / (1. + gamma);
}
 
template <typename T0, typename T1, typename T2>
void l2_norm(Eigen::DenseBase<T0> &out, Eigen::DenseBase<T2> const &gamma,
             Eigen::DenseBase<T1> const &x) {
  out = x.derived().array() * 1. / (1. + gamma.derived().array()).array();
}
 
template <typename T0, typename T1>
void tv_norm(Eigen::DenseBase<T0> &out, typename real_type<typename T0::Scalar>::type gamma,
             Eigen::DenseBase<T1> const &x) {
  typename Eigen::Index const size = x.size() / 2;
  typename T1::PlainObject const &u = x.segment(0, size);
  typename T1::PlainObject const &v = x.segment(size, size);
  out = T0::Zero(size * 2);
  for (typename Eigen::Index i(0); i < size; i++) {
    const t_real norm = std::sqrt(std::abs(u(i) * u(i) + v(i) * v(i)));
    if (norm > gamma) {
      out(i) = (1 - gamma / norm) * u(i);
      out(size + i) = (1 - gamma / norm) * v(i);
    } else {
      out(i) = 0;
      out(size + i) = 0;
    }
  }
}
template <typename T0, typename T1, typename T2>
void tv_norm(Eigen::DenseBase<T0> &out, Eigen::DenseBase<T2> const &gamma,
             Eigen::DenseBase<T1> const &x) {
  typename Eigen::Index const size = x.size() / 2;
  typename T1::PlainObject const &u = x.segment(0, size);
  typename T1::PlainObject const &v = x.segment(size, size);
  out = T0::Zero(size * 2);
  for (typename Eigen::Index i(0); i < size; i++) {
    const t_real norm = std::sqrt(std::abs(u(i) * u(i) + v(i) * v(i)));
    if (norm > gamma(i)) {
      out(i) = (1 - gamma(i) / norm) * u(i);
      out(size + i) = (1 - gamma(i) / norm) * v(i);
    } else {
      out(i) = 0;
      out(size + i) = 0;
    }
  }
}
 
template <typename T0, typename T1>
void id(Eigen::DenseBase<T0> &out, typename real_type<typename T0::Scalar>::type gamma,
        Eigen::DenseBase<T1> const &x) {
  out = x;
}
 
template <typename T>
auto l1_norm(typename real_type<typename T::Scalar>::type gamma, Eigen::DenseBase<T> const &x)
    -> decltype(sopt::soft_threshhold(x, gamma)) {
  return sopt::soft_threshhold(x, gamma);
}
 
template <typename T>
void positive_quadrant(Vector<T> &out, typename real_type<T>::type, Vector<T> const &x) {
  out = sopt::positive_quadrant(x);
};
 
template <typename T>
class L2Norm {
 public:
  using Real = typename real_type<T>::type;
  L2Norm() {}
 
  void operator()(Vector<T> &out, const Real gamma, Vector<T> const &x) const {
    proximal::l2_norm(out, gamma, x);
  }
  template <typename T0>
  EnveloppeExpression<L2Norm, T0> operator()(Real const &, Eigen::MatrixBase<T0> const &x) const {
    return {*this, x};
  }
 
  template <typename T0>
  EnveloppeExpression<L2Norm, T0> operator()(Eigen::MatrixBase<T0> const &x) const {
    return {*this, x};
  }
};
 
template <typename T>
class L2Ball {
 public:
  using Real = typename real_type<T>::type;
#ifdef SOPT_MPI
  L2Ball(Real epsilon, mpi::Communicator const &comm = mpi::Communicator())
      : epsilon_(epsilon), comm_(comm){};
#else
  L2Ball(Real epsilon) : epsilon_(epsilon) {}
#endif
 
  void operator()(Vector<T> &out, typename real_type<T>::type, Vector<T> const &x) const {
    return operator()(out, x);
  }
  void operator()(Vector<T> &out, Vector<T> const &x) const {
#ifdef SOPT_MPI
    auto const norm = mpi::l2_norm(x, comm_);
#else
    auto const norm = x.stableNorm();
#endif
    if (norm > epsilon())
      out = x * (epsilon() / norm);
    else
      out = x;
  }
  template <typename T0>
  EnveloppeExpression<L2Ball, T0> operator()(Real const &, Eigen::MatrixBase<T0> const &x) const {
    return {*this, x};
  }
 
  template <typename T0>
  EnveloppeExpression<L2Ball, T0> operator()(Eigen::MatrixBase<T0> const &x) const {
    return {*this, x};
  }
 
  Real epsilon() const { return epsilon_; }
  L2Ball<T> &epsilon(Real eps) {
    epsilon_ = eps;
    return *this;
  }
 
#ifdef SOPT_MPI
  mpi::Communicator const &communicator() const { return comm_; }
  L2Ball &communicator(mpi::Communicator const &comm) {
    comm_ = comm;
    return *this;
  }
#endif
 
 private:
  Real epsilon_;
#ifdef SOPT_MPI
  mpi::Communicator comm_;
#endif
};
 
template <typename T>
class WeightedL2Ball : public L2Ball<T> {
 public:
  using Real = typename L2Ball<T>::Real;
  using t_Vector = Vector<Real>;
#ifdef SOPT_MPI
  template <typename T0>
  WeightedL2Ball(Real epsilon, Eigen::DenseBase<T0> const &w,
                 mpi::Communicator const &comm = mpi::Communicator())
      : L2Ball<T>(epsilon, comm), weights_(w) {}
  WeightedL2Ball(Real epsilon, mpi::Communicator const &comm = mpi::Communicator())
      : WeightedL2Ball(epsilon, t_Vector::Ones(1), comm) {}
  mpi::Communicator communicator() const { return L2Ball<T>::communicator(); }
  WeightedL2Ball<T> &communicator(mpi::Communicator const &c) {
    L2Ball<T>::communicator(c);
    return *this;
  }
#else
  template <typename T0>
  WeightedL2Ball(Real epsilon, Eigen::DenseBase<T0> const &w) : L2Ball<T>(epsilon), weights_(w) {}
  WeightedL2Ball(Real epsilon) : WeightedL2Ball(epsilon, t_Vector::Ones(1)) {}
#endif
 
  void operator()(Vector<T> &out, typename real_type<T>::type, Vector<T> const &x) const {
    return operator()(out, x);
  }
  void operator()(Vector<T> &out, Vector<T> const &x) const {
#ifdef SOPT_MPI
    auto const norm = mpi::l2_norm(x.array(), weights().array(), communicator());
#else
    auto const norm = weights().size() == 1 ? x.stableNorm() * std::abs(weights()(0))
                                            : (x.array() * weights().array()).matrix().stableNorm();
#endif
    if (norm > epsilon())
      out = x * (epsilon() / norm);
    else
      out = x;
  }
  template <typename T0>
  EnveloppeExpression<WeightedL2Ball, T0> operator()(Real const &,
                                                     Eigen::MatrixBase<T0> const &x) const {
    return {*this, x};
  }
  template <typename T0>
  EnveloppeExpression<WeightedL2Ball, T0> operator()(Eigen::MatrixBase<T0> const &x) const {
    return {*this, x};
  }
 
  t_Vector const &weights() const { return weights_; }
  template <typename T0>
  WeightedL2Ball<T> &weights(Eigen::MatrixBase<T0> const &w) {
    if ((w.array() < 0e0).any()) SOPT_THROW("Weights cannot be negative");
    if (w.stableNorm() < 1e-12) SOPT_THROW("Weights cannot be null");
    weights_ = w;
    return *this;
  }
  Real epsilon() const { return L2Ball<T>::epsilon(); }
  WeightedL2Ball<T> &epsilon(Real const &eps) {
    L2Ball<T>::epsilon(eps);
    return *this;
  }
 
 private:
  t_Vector weights_;
};
 
template <typename FUNCTION, typename VECTOR>
class Translation {
 public:
  template <typename T_VECTOR>
  Translation(FUNCTION const &func, T_VECTOR const &trans) : func(func), trans(trans) {}
  template <typename OUTPUT, typename T0>
  typename std::enable_if<std::is_reference<OUTPUT>::value, void>::type operator()(
      OUTPUT out, typename real_type<typename T0::Scalar>::type const &t,
      Eigen::MatrixBase<T0> const &x) const {
    func(out, t, x + trans);
    out -= trans;
  }
  template <typename T0>
  void operator()(Vector<typename T0::Scalar> &out,
                  typename real_type<typename T0::Scalar>::type const &t,
                  Eigen::MatrixBase<T0> const &x) const {
    func(out, t, x + trans);
    out -= trans;
  }
  template <typename T0>
  ProximalExpression<Translation<FUNCTION, VECTOR> const &, T0> operator()(
      typename T0::Scalar const &t, Eigen::MatrixBase<T0> const &x) const {
    return {*this, t, x};
  }
 
 private:
  FUNCTION const func;
  VECTOR const trans;
};
 
template <typename FUNCTION, typename VECTOR>
Translation<FUNCTION, VECTOR> translate(FUNCTION const &func, VECTOR const &translation) {
  return Translation<FUNCTION, VECTOR>(func, translation);
}
} // namespace sopt::proximal
 
#endif