Merge pull request #84 from mrava87/main

Variable sigma in L1

Author: mrava87

pyproximal/optimization/palm.py

@@ -77,28 +77,33 @@ def PALM(H, proxf, proxg, x0, y0, gammaf=1., gammag=1.,
               'Proximal operator (g): %s\n'
               'gammaf = %10e\tgammaf = %10e\tniter = %d\n' %
               (type(H), type(proxf), type(proxg), gammaf, gammag, niter))
-        head = '   Itn      x[0]       y[0]        f         g         H'
+        head = '   Itn      x[0]       y[0]        f         g         H         ck         dk'
         print(head)
 
     x, y = x0.copy(), y0.copy()
     for iiter in range(niter):
         ck = gammaf * H.ly(y)
         x = x - (1 / ck) * H.gradx(x.ravel())
-        x = proxf.prox(x, ck)
+        if proxf is not None:
+            x = proxf.prox(x, ck)
         H.updatex(x.copy())
         dk = gammag * H.lx(x)
         y = y - (1 / dk) * H.grady(y.ravel())
-        y = proxg.prox(y, dk)
+        if proxg is not None:
+            y = proxg.prox(y, dk)
         H.updatey(y.copy())
 
         # run callback
         if callback is not None:
             callback(x, y)
 
         if show:
+            pf = proxf(x) if proxf is not None else 0.
+            pg = proxg(y) if proxg is not None else 0.
             if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0:
-                msg = '%6g  %5.5e  %5.2e  %5.2e  %5.2e  %5.2e' % \
-                      (iiter + 1, x[0], y[0], proxf(x), proxg(y), H(x, y))
+                msg = '%6g  %5.5e  %5.2e  %5.2e  %5.2e  %5.2e  %5.2e  %5.2e' % \
+                      (iiter + 1, x[0], y[0], pf if pf is not None else 0.,
+                       pg if pg is not None else 0., H(x, y), ck, dk)
                 print(msg)
     if show:
         print('\nTotal time (s) = %.2f' % (time.time() - tstart))

pyproximal/proximal/L1.py

@@ -33,6 +33,13 @@ def _softthreshold(x, thresh):
     return x1
 
 
+def _current_sigma(sigma, count):
+    if not callable(sigma):
+        return sigma
+    else:
+        return sigma(count)
+
+
 class L1(ProxOperator):
     r"""L1 norm proximal operator.
 
@@ -41,8 +48,12 @@ class L1(ProxOperator):
 
     Parameters
     ----------
-    sigma : :obj:`int`, optional
-        Multiplicative coefficient of L1 norm
+    sigma : :obj:`float` or :obj:`list` or :obj:`np.ndarray` or :obj:`func`, optional
+        Multiplicative coefficient of L1 norm. This can be a constant number, a list
+        of values (for multidimensional inputs, acting on the second dimension) or
+        a function that is called passing a counter which keeps track of how many
+        times the ``prox`` method has been invoked before and returns a scalar (or a list of)
+        ``sigma`` to be used.
     g : :obj:`np.ndarray`, optional
         Vector to be subtracted
 
@@ -85,18 +96,33 @@ def __init__(self, sigma=1., g=None):
         self.sigma = sigma
         self.g = g
         self.gdual = 0 if g is None else g
-        self.box = BoxProj(-sigma, sigma)
+        if not callable(sigma):
+            self.box = BoxProj(-sigma, sigma)
+        else:
+            self.box = BoxProj(-sigma(0), sigma(0))
+        self.count = 0
 
     def __call__(self, x):
-        return self.sigma * np.sum(np.abs(x))
-
+        sigma = _current_sigma(self.sigma, self.count)
+        return sigma * np.sum(np.abs(x))
+
+    def _increment_count(func):
+        """Increment counter
+        """
+        def wrapped(self, *args, **kwargs):
+            self.count += 1
+            return func(self, *args, **kwargs)
+        return wrapped
+
+    @_increment_count
     @_check_tau
     def prox(self, x, tau):
+        sigma = _current_sigma(self.sigma, self.count)
         if self.g is None:
-            x = _softthreshold(x, tau * self.sigma)
+            x = _softthreshold(x, tau * sigma)
         else:
             # use precomposition property
-            x = _softthreshold(x - self.g, tau * self.sigma) + self.g
+            x = _softthreshold(x - self.g, tau * sigma) + self.g
         return x
 
     @_check_tau

pyproximal/proximal/L2.py

@@ -97,7 +97,7 @@ def __init__(self, Op=None, b=None, q=None, sigma=1., alpha=1.,
         self.count = 0
 
         # create data term
-        if self.Op is not None:
+        if self.Op is not None and self.b is not None:
             self.OpTb = self.sigma * self.Op.H @ self.b
             # create A.T A upfront for explicit operators
             if self.Op.explicit: