astro-informatics · tkoskela · Jan 9, 2023 · Dec 8, 2022 · Dec 8, 2022 · Dec 9, 2022
diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml
@@ -77,7 +77,8 @@ jobs:
     - name: Install Dependencies on MacOS
       if: ${{ contains(matrix.os, 'macos') }}
       run: |
-        brew update
+        # Skip brew update until https://github.com/actions/setup-python/issues/577 is fixed
+        # brew update
         brew install libtiff open-mpi libyaml ccache
         echo "CMAKE_PREFIX_PATH=/usr/local/opt/libomp" >> $GITHUB_ENV
 
@@ -87,14 +88,14 @@ jobs:
       uses: actions/setup-python@v4
       with:
         python-version: '3.10'
-        
+
     - name: Install Conan
       id: conan
       uses: turtlebrowser/get-conan@main
 
     - name: Conan version
       run: echo "${{ steps.conan.outputs.version }}"
-        
+
     - name: Prepare ccache timestamp
       id: ccache_cache_timestamp
       run: echo "{date_and_time}={$(date +'%Y-%m-%d-%H;%M;%S')}" >> $GITHUB_OUTPUT

diff --git a/cpp/examples/forward_backward/inpainting.cc b/cpp/examples/forward_backward/inpainting.cc
@@ -106,11 +106,13 @@ int main(int argc, char const **argv) {
   auto gp = std::make_shared<sopt::algorithm::L1GProximal<Scalar>>(false);
   gp->l1_proximal_tolerance(1e-4)
     .l1_proximal_nu(1)
+    .l1_proximal_beta(fb.beta())
+    .l1_proximal_Phi(fb.Phi())
     .l1_proximal_itermax(50)
     .l1_proximal_positivity_constraint(true)
     .l1_proximal_real_constraint(true)
     .Psi(psi);
-  
+
   // Once the properties are set, inject it into the ImagingForwardBackward object
   fb.g_proximal(gp);
 

diff --git a/cpp/examples/forward_backward/inpainting_credible_interval.cc b/cpp/examples/forward_backward/inpainting_credible_interval.cc
@@ -108,14 +108,16 @@ int main(int argc, char const **argv) {
   auto gp = std::make_shared<sopt::algorithm::L1GProximal<Scalar>>(false);
   gp->l1_proximal_tolerance(1e-4)
     .l1_proximal_nu(1)
+    .l1_proximal_beta(fb.beta())
+    .l1_proximal_Phi(fb.Phi())
     .l1_proximal_itermax(50)
     .l1_proximal_positivity_constraint(true)
     .l1_proximal_real_constraint(true)
     .Psi(psi);
-  
+
   // Once the properties are set, inject it into the ImagingForwardBackward object
   fb.g_proximal(gp);
-  
+
   SOPT_HIGH_LOG("Starting Forward Backward");
   // Alternatively, forward-backward can be called with a tuple (x, residual) as argument
   // Here, we default to (Φ^Ty/ν, ΦΦ^Ty/ν - y)

diff --git a/cpp/examples/forward_backward/inpainting_joint_map.cc b/cpp/examples/forward_backward/inpainting_joint_map.cc
@@ -107,11 +107,13 @@ int main(int argc, char const **argv) {
   auto gp = std::make_shared<sopt::algorithm::L1GProximal<Scalar>>(false);
   gp->l1_proximal_tolerance(1e-4)
     .l1_proximal_nu(1)
+    .l1_proximal_beta(fb->beta())
+    .l1_proximal_Phi(fb->Phi())
     .l1_proximal_itermax(50)
     .l1_proximal_positivity_constraint(true)
     .l1_proximal_real_constraint(true)
     .Psi(psi);
-  
+
   // Once the properties are set, inject it into the ImagingForwardBackward object
   fb->g_proximal(gp);
 

diff --git a/cpp/sopt/forward_backward.h b/cpp/sopt/forward_backward.h
@@ -231,7 +231,7 @@ void ForwardBackward<SCALAR>::iteration_step(t_Vector &out, t_Vector &residual,
                                              t_Vector &z, const t_real lambda) const {
   p = out;
   f_gradient(z, residual);
-  g_proximal(out, gamma() * beta(), out - beta() / nu() * (Phi().adjoint() * z));
+  g_proximal(out, gamma(), z);
   p = out + lambda * (out - p);
   residual = (Phi() * p) - target();
 }

diff --git a/cpp/sopt/l1_g_proximal.h b/cpp/sopt/l1_g_proximal.h
@@ -55,7 +55,9 @@ class L1GProximal : public GProximal<SCALAR> {
   // Return g_proximal as a lambda function. Used in operator() in base class.
   t_Proximal proximal_function() const override {
     return [this](t_Vector &out, Real gamma, t_Vector const &x) {
-	     this -> l1_proximal(out, gamma, x);
+	     this -> l1_proximal(out,
+				 gamma * l1_proximal_.beta(),
+				 out - l1_proximal_.beta() / l1_proximal_.nu() * (l1_proximal_.Phi().adjoint() * x));
 	   };
   }
 
@@ -85,10 +87,10 @@ class L1GProximal : public GProximal<SCALAR> {
 // ~~~
 #define SOPT_MACRO(VAR, TYPE)						 \
   /** \brief Getter, forwards to l1_proximal **/                         \
-  TYPE const &l1_proximal_##VAR() const { return l1_proximal().VAR(); }  \
+  TYPE const &l1_proximal_##VAR() const { return l1_proximal_.VAR(); }  \
   /** \brief Setter, forwards to l1_proximal **/                         \
   L1GProximal<SCALAR> &l1_proximal_##VAR(TYPE const ARG) {               \
-    l1_proximal().VAR(ARG);				                 \
+    l1_proximal_.VAR(ARG);				                 \
     return *this;                                                        \
   }
 
@@ -97,7 +99,9 @@ class L1GProximal : public GProximal<SCALAR> {
   SOPT_MACRO(positivity_constraint, bool);
   SOPT_MACRO(real_constraint, bool);
   SOPT_MACRO(nu, Real);
+  SOPT_MACRO(beta, Real);
   SOPT_MACRO(weights, t_Vector);
+  SOPT_MACRO(Phi, t_LinearTransform);
 #undef SOPT_MACRO
 
   //! Analysis operator Ψ

diff --git a/cpp/sopt/l1_proximal.h b/cpp/sopt/l1_proximal.h
@@ -40,12 +40,16 @@ class L1TightFrame {
                mpi::Communicator const &adjoint_comm = mpi::Communicator())
       : Psi_(linear_transform_identity<Scalar>()),
         nu_(1e0),
+	Phi_(linear_transform_identity<Scalar>()),
+	beta_(1e0),
         direct_space_comm_(direct_comm),
         adjoint_space_comm_(adjoint_comm),
         weights_(Vector<Real>::Ones(1)) {}
 #else
   L1TightFrame()
-      : Psi_(linear_transform_identity<Scalar>()), nu_(1e0), weights_(Vector<Real>::Ones(1)) {}
+    : Psi_(linear_transform_identity<Scalar>()), nu_(1e0),
+      Phi_(linear_transform_identity<Scalar>()), beta_(1e0),
+      weights_(Vector<Real>::Ones(1))  {}
 #endif
 
 #define SOPT_MACRO(NAME, TYPE)                   \
@@ -62,7 +66,11 @@ class L1TightFrame {
   //! Linear transform applied to input prior to L1 norm
   SOPT_MACRO(Psi, LinearTransform<Vector<Scalar>>);
   //! Bound on the squared norm of the operator Ψ
-  SOPT_MACRO(nu, Real);
+  SOPT_MACRO(nu, Real)
+  //! Measurement operator, used by g_proximal
+  SOPT_MACRO(Phi, LinearTransform<Vector<Scalar>>);
+  //! β parameter, used by g_proximal
+  SOPT_MACRO(beta, Real);
 #ifdef SOPT_MPI
   //! Communicator for summing in direct space (input when applying Psi)
   SOPT_MACRO(direct_space_comm, mpi::Communicator);
@@ -307,6 +315,24 @@ class L1 : protected L1TightFrame<SCALAR> {
     return *this;
   }
 
+  //! Bounds on the squared norm of the operator Ψ
+  Real beta() const { return L1TightFrame<Scalar>::beta(); }
+  //! Sets the bound on the squared norm of the operator Ψ
+  L1<Scalar> &beta(Real const &beta) {
+    L1TightFrame<SCALAR>::beta(beta);
+    return *this;
+  }
+
+  //! Linear transform applied to input prior to L1 norm
+  LinearTransform<Vector<Scalar>> const &Phi() const { return L1TightFrame<Scalar>::Phi(); }
+  //! Set Ψ and Ψ^† using a matrix
+  template <class... ARGS>
+  typename std::enable_if<sizeof...(ARGS) >= 1, L1<Scalar> &>::type Phi(ARGS &&... args) {
+    L1TightFrame<Scalar>::Phi(std::forward<ARGS>(args)...);
+    return *this;
+  }
+
+
   //! \brief Special case if Ψ ia a tight frame.
   //! \see L1TightFrame
   template <class... T>

diff --git a/cpp/sopt/l2_forward_backward.h b/cpp/sopt/l2_forward_backward.h
@@ -259,9 +259,11 @@ typename L2ForwardBackward<SCALAR>::Diagnostic L2ForwardBackward<SCALAR>::operat
   Diagnostic result;
   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->l2_proximal_weights() * gamma, x);
+      this->l2_proximal_weighted()(out, this->l2_proximal_weights() * gamma,
+				   out - beta() / nu() * (Phi().adjoint() * x));
     else
-      this->l2_proximal()(out, this->l2_proximal_weights()(0) * gamma, x);
+      this->l2_proximal()(out, this->l2_proximal_weights()(0) * gamma,
+			  out - beta() / nu() * (Phi().adjoint() * x));
   };
   const Real sigma_factor = sigma() * sigma();
   auto const f_gradient = [this, sigma_factor](t_Vector &out, t_Vector const &x) {

diff --git a/cpp/tests/forward_backward.cc b/cpp/tests/forward_backward.cc
@@ -25,13 +25,19 @@ sopt::t_int random_integer(sopt::t_int min, sopt::t_int max) {
 typedef sopt::t_real Scalar;
 typedef sopt::Vector<Scalar> t_Vector;
 typedef sopt::t_real t_real;
+typedef sopt::LinearTransform<t_Vector> t_LinearTransform;
 auto const N = 5;
 
 TEST_CASE("Forward Backward with ||x - x0||_2^2 function", "[fb]") {
   using namespace sopt;
   t_Vector const target0 = t_Vector::Random(N);
-  auto const g0 = [](t_Vector &out, const t_real gamma, const t_Vector &x) {
-    proximal::id(out, gamma, x);
+  t_real const beta = 0.2;
+  t_real const nu = 1.0;
+  t_real const gamma = 0.1;
+  int const itermax = 300;
+  t_LinearTransform const Phi = linear_transform_identity<Scalar>();
+  auto const g0 = [=](t_Vector &out, const t_real gamma, const t_Vector &x) {
+		    proximal::id(out, gamma * beta, out - beta / nu * (Phi.adjoint() * x));
   };
   auto const grad = [](t_Vector &out, const t_Vector &x) { out = x; };
   const t_Vector x_guess = t_Vector::Random(target0.size());
@@ -43,10 +49,11 @@ TEST_CASE("Forward Backward with ||x - x0||_2^2 function", "[fb]") {
   CAPTURE(x_guess);
   CAPTURE(res);
   auto const fb = algorithm::ForwardBackward<Scalar>(grad, g0, target0)
-                      .itermax(300)
-                      .gamma(0.1)
-                      .beta(0.2)
+                      .itermax(itermax)
+                      .gamma(gamma)
+                      .beta(beta)
                       .is_converged(convergence);
+
   auto const result = fb(std::make_tuple(x_guess, res));
   CAPTURE(result.niters);
   CAPTURE(result.x);
@@ -86,6 +93,8 @@ TEST_CASE("Check type returned on setting variables") {
   auto gp = std::make_shared<sopt::algorithm::L1GProximal<Scalar>>(false);
   CHECK(is_l1_g_proximal_ref<decltype(gp->l1_proximal_tolerance(1e-2))>::value);
   CHECK(is_l1_g_proximal_ref<decltype(gp->l1_proximal_nu(1))>::value);
+  CHECK(is_l1_g_proximal_ref<decltype(gp->l1_proximal_beta(1))>::value);
+  CHECK(is_l1_g_proximal_ref<decltype(gp->l1_proximal_Phi(linear_transform_identity<double>()))>::value);
   CHECK(is_l1_g_proximal_ref<decltype(gp->l1_proximal_itermax(50))>::value);
   CHECK(is_l1_g_proximal_ref<decltype(gp->l1_proximal_positivity_constraint(true))>::value);
   CHECK(is_l1_g_proximal_ref<decltype(gp->l1_proximal_real_constraint(true))>::value);

diff --git a/cpp/tests/inpainting.cc b/cpp/tests/inpainting.cc
@@ -70,11 +70,13 @@ TEST_CASE("Inpainting"){
   auto gp = std::make_shared<sopt::algorithm::L1GProximal<Scalar>>(false);
   gp->l1_proximal_tolerance(1e-4)
     .l1_proximal_nu(1)
+    .l1_proximal_beta(fb.beta())
+    .l1_proximal_Phi(fb.Phi())
     .l1_proximal_itermax(50)
     .l1_proximal_positivity_constraint(true)
     .l1_proximal_real_constraint(true)
     .Psi(psi);
-  
+
   // Once the properties are set, inject it into the ImagingForwardBackward object
   fb.g_proximal(gp);