diff --git a/docs/source/api/index.rst b/docs/source/api/index.rst index 8be754b..b8b92be 100755 --- a/docs/source/api/index.rst +++ b/docs/source/api/index.rst @@ -64,6 +64,7 @@ Convex EuclideanBall Huber Intersection + L0 L0Ball L1 L1Ball @@ -136,6 +137,7 @@ Primal AcceleratedProximalGradient ADMM ADMML2 + HQS LinearizedADMM ProximalGradient ProximalPoint diff --git a/pyproximal/optimization/__init__.py b/pyproximal/optimization/__init__.py index 5b8f5ab..c25ac22 100644 --- a/pyproximal/optimization/__init__.py +++ b/pyproximal/optimization/__init__.py @@ -11,8 +11,9 @@ ProximalPoint Proximal point algorithm (or proximal min.) ProximalGradient Proximal gradient algorithm AcceleratedProximalGradient Accelerated Proximal gradient algorithm + HQS Half Quadrating Splitting ADMM Alternating Direction Method of Multipliers - ADMML2 ADMM with L2 misfit term + ADMML2 ADMM with L2 misfit term LinearizedADMM Linearized ADMM TwIST Two-step Iterative Shrinkage/Threshold PlugAndPlay Plug-and-Play Prior with ADMM diff --git a/pyproximal/optimization/primal.py b/pyproximal/optimization/primal.py index 4021ed8..cd4e8ab 100644 --- a/pyproximal/optimization/primal.py +++ b/pyproximal/optimization/primal.py @@ -372,7 +372,115 @@ def AcceleratedProximalGradient(proxf, proxg, x0, tau=None, beta=0.5, return x -def ADMM(proxf, proxg, x0, tau, niter=10, callback=None, show=False): +def HQS(proxf, proxg, x0, tau, niter=10, gfirst=True, + callback=None, callbackz=False, show=False): + r"""Half Quadratic splitting + + Solves the following minimization problem using Half Quadratic splitting + algorithm: + + .. math:: + + \mathbf{x},\mathbf{z} = \argmin_{\mathbf{x},\mathbf{z}} + f(\mathbf{x}) + g(\mathbf{z}) \\ + s.t. \; \mathbf{x}=\mathbf{z} + + where :math:`f(\mathbf{x})` and :math:`g(\mathbf{z})` are any convex + function that has a known proximal operator. + + Parameters + ---------- + proxf : :obj:`pyproximal.ProxOperator` + Proximal operator of f function + proxg : :obj:`pyproximal.ProxOperator` + Proximal operator of g function + x0 : :obj:`numpy.ndarray` + Initial vector + tau : :obj:`float`, optional + Positive scalar weight, which should satisfy the following condition + to guarantees convergence: :math:`\tau \in (0, 1/L]` where ``L`` is + the Lipschitz constant of :math:`\nabla f`. + niter : :obj:`int`, optional + Number of iterations of iterative scheme + gfirst : :obj:`bool`, optional + Apply Proximal of operator ``g`` first (``True``) or Proximal of + operator ``f`` first (``False``) + callback : :obj:`callable`, optional + Function with signature (``callback(x)``) to call after each iteration + where ``x`` is the current model vector + callbackz : :obj:`bool`, optional + Modify callback signature to (``callback(x, z)``) when ``callbackz=True`` + show : :obj:`bool`, optional + Display iterations log + + Returns + ------- + x : :obj:`numpy.ndarray` + Inverted model + z : :obj:`numpy.ndarray` + Inverted second model + + Notes + ----- + The HQS algorithm can be expressed by the following recursion [1]_: + + .. math:: + + \mathbf{z}^{k+1} = \prox_{\tau g}(\mathbf{x}^{k}) + \mathbf{x}^{k+1} = \prox_{\tau f}(\mathbf{z}^{k+1})\\ + + Note that ``x`` and ``z`` converge to each other, however if iterations are + stopped too early ``x`` is guaranteed to belong to the domain of ``f`` + while ``z`` is guaranteed to belong to the domain of ``g``. Depending on + the problem either of the two may be the best solution. + + .. [1] D., Geman, and C., Yang, "Nonlinear image recovery with halfquadratic + regularization", IEEE Transactions on Image Processing, + 4, 7, pp. 932-946, 1995. + + """ + if show: + tstart = time.time() + print('HQS\n' + '---------------------------------------------------------\n' + 'Proximal operator (f): %s\n' + 'Proximal operator (g): %s\n' + 'tau = %10e\tniter = %d\n' % (type(proxf), type(proxg), + tau, niter)) + head = ' Itn x[0] f g J = f + g' + print(head) + + x = x0.copy() + z = np.zeros_like(x) + for iiter in range(niter): + if gfirst: + z = proxg.prox(x, tau) + x = proxf.prox(z, tau) + else: + x = proxf.prox(z, tau) + z = proxg.prox(x, tau) + + # run callback + if callback is not None: + if callbackz: + callback(x, z) + else: + callback(x) + + if show: + if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0: + pf, pg = proxf(x), proxg(x) + msg = '%6g %12.5e %10.3e %10.3e %10.3e' % \ + (iiter + 1, x[0], pf, pg, pf + pg) + print(msg) + if show: + print('\nTotal time (s) = %.2f' % (time.time() - tstart)) + print('---------------------------------------------------------\n') + return x, z + + +def ADMM(proxf, proxg, x0, tau, niter=10, gfirst=False, + callback=None, callbackz=False, show=False): r"""Alternating Direction Method of Multipliers Solves the following minimization problem using Alternating Direction @@ -417,9 +525,14 @@ def ADMM(proxf, proxg, x0, tau, niter=10, callback=None, show=False): the Lipschitz constant of :math:`\nabla f`. niter : :obj:`int`, optional Number of iterations of iterative scheme + gfirst : :obj:`bool`, optional + Apply Proximal of operator ``g`` first (``True``) or Proximal of + operator ``f`` first (``False``) callback : :obj:`callable`, optional Function with signature (``callback(x)``) to call after each iteration where ``x`` is the current model vector + callbackz : :obj:`bool`, optional + Modify callback signature to (``callback(x, z)``) when ``callbackz=True`` show : :obj:`bool`, optional Display iterations log @@ -445,7 +558,7 @@ def ADMM(proxf, proxg, x0, tau, niter=10, callback=None, show=False): \mathbf{z}^{k+1} = \prox_{\tau g}(\mathbf{x}^{k+1} + \mathbf{u}^{k})\\ \mathbf{u}^{k+1} = \mathbf{u}^{k} + \mathbf{x}^{k+1} - \mathbf{z}^{k+1} - Note that ``x`` and ``z`` converge to each other, but if iterations are + Note that ``x`` and ``z`` converge to each other, however if iterations are stopped too early ``x`` is guaranteed to belong to the domain of ``f`` while ``z`` is guaranteed to belong to the domain of ``g``. Depending on the problem either of the two may be the best solution. @@ -465,14 +578,20 @@ def ADMM(proxf, proxg, x0, tau, niter=10, callback=None, show=False): x = x0.copy() u = z = np.zeros_like(x) for iiter in range(niter): - x = proxf.prox(z - u, tau) - z = proxg.prox(x + u, tau) + if gfirst: + z = proxg.prox(x + u, tau) + x = proxf.prox(z - u, tau) + else: + x = proxf.prox(z - u, tau) + z = proxg.prox(x + u, tau) u = u + x - z # run callback if callback is not None: - callback(x) - + if callbackz: + callback(x, z) + else: + callback(x) if show: if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0: pf, pg = proxf(x), proxg(x) diff --git a/pyproximal/optimization/primaldual.py b/pyproximal/optimization/primaldual.py index db58635..99e3dca 100644 --- a/pyproximal/optimization/primaldual.py +++ b/pyproximal/optimization/primaldual.py @@ -3,7 +3,7 @@ def PrimalDual(proxf, proxg, A, x0, tau, mu, z=None, theta=1., niter=10, - gfirst=True, callback=None, show=False): + gfirst=True, callback=None, callbacky=False, show=False): r"""Primal-dual algorithm Solves the following (possibly) nonlinear minimization problem using @@ -39,10 +39,12 @@ def PrimalDual(proxf, proxg, A, x0, tau, mu, z=None, theta=1., niter=10, Linear operator of g x0 : :obj:`numpy.ndarray` Initial vector - tau : :obj:`float` - Stepsize of subgradient of :math:`f` - mu : :obj:`float` - Stepsize of subgradient of :math:`g^*` + tau : :obj:`float` or :obj:`np.ndarray` + Stepsize of subgradient of :math:`f`. This can be constant + or function of iterations (in the latter cases provided as np.ndarray) + mu : :obj:`float` or :obj:`np.ndarray` + Stepsize of subgradient of :math:`g^*`. This can be constant + or function of iterations (in the latter cases provided as np.ndarray) z : :obj:`numpy.ndarray`, optional Additional vector theta : :obj:`float` @@ -58,6 +60,8 @@ def PrimalDual(proxf, proxg, A, x0, tau, mu, z=None, theta=1., niter=10, callback : :obj:`callable`, optional Function with signature (``callback(x)``) to call after each iteration where ``x`` is the current model vector + callbacky : :obj:`bool`, optional + Modify callback signature to (``callback(x, y)``) when ``callbacky=True`` show : :obj:`bool`, optional Display iterations log @@ -93,11 +97,20 @@ def PrimalDual(proxf, proxg, A, x0, tau, mu, z=None, theta=1., niter=10, \mathbf{y}^{k+1} = \prox_{\mu g^*}(\mathbf{y}^{k} + \mu \mathbf{A}\bar{\mathbf{x}}^{k+1}) - .. [1] A., Chambolle, and T., Pock, "A first-order primal-dual algorithm for + .. [1] 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, 8pp. 120–145. 2011. + Imaging and Vision, 40, 8pp. 120-145. 2011. """ + # check if tau and mu are scalars or arrays + fixedtau = fixedmu = False + if isinstance(tau, (int, float)): + tau = tau * np.ones(niter) + fixedtau = True + if isinstance(mu, (int, float)): + mu = mu * np.ones(niter) + fixedmu = True + if show: tstart = time.time() print('Primal-dual: min_x f(Ax) + x^T z + g(x)\n' @@ -106,9 +119,10 @@ def PrimalDual(proxf, proxg, A, x0, tau, mu, z=None, theta=1., niter=10, 'Proximal operator (g): %s\n' 'Linear operator (A): %s\n' 'Additional vector (z): %s\n' - 'tau = %10e\tmu = %10e\ntheta = %.2f\t\tniter = %d\n' % + 'tau = %s\t\tmu = %s\ntheta = %.2f\t\tniter = %d\n' % (type(proxf), type(proxg), type(A), - None if z is None else 'vector', tau, mu, theta, niter)) + None if z is None else 'vector', str(tau[0]) if fixedtau else 'Variable', + str(mu[0]) if fixedmu else 'Variable', theta, niter)) head = ' Itn x[0] f g z^x J = f + g + z^x' print(head) @@ -119,24 +133,26 @@ def PrimalDual(proxf, proxg, A, x0, tau, mu, z=None, theta=1., niter=10, for iiter in range(niter): xold = x.copy() if gfirst: - y = proxg.proxdual(y + mu * A.matvec(xhat), mu) + y = proxg.proxdual(y + mu[iiter] * A.matvec(xhat), mu[iiter]) ATy = A.rmatvec(y) if z is not None: ATy += z - x = proxf.prox(x - tau * ATy, tau) + x = proxf.prox(x - tau[iiter] * ATy, tau[iiter]) xhat = x + theta * (x - xold) else: ATy = A.rmatvec(y) if z is not None: ATy += z - x = proxf.prox(x - tau * ATy, tau) + x = proxf.prox(x - tau[iiter] * ATy, tau[iiter]) xhat = x + theta * (x - xold) - y = proxg.proxdual(y + mu * A.matvec(xhat), mu) + y = proxg.proxdual(y + mu[iiter] * A.matvec(xhat), mu[iiter]) # run callback if callback is not None: - callback(x) - + if callbacky: + callback(x, y) + else: + callback(x) if show: if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0: pf, pg = proxf(x), proxg(A.matvec(x)) diff --git a/pyproximal/projection/L1.py b/pyproximal/projection/L1.py index d4a1ed8..cfd4b3c 100644 --- a/pyproximal/projection/L1.py +++ b/pyproximal/projection/L1.py @@ -40,4 +40,8 @@ def __init__(self, n, radius, maxiter=100, xtol=1e-5): self.simplex = SimplexProj(n, radius, maxiter, xtol) def __call__(self, x): - return np.sign(x) * self.simplex(np.abs(x)) \ No newline at end of file + if np.iscomplexobj(x): + return np.exp(1j * np.angle(x)) * self.simplex(np.abs(x)) + else: + return np.sign(x) * self.simplex(np.abs(x)) + diff --git a/pyproximal/proximal/L0.py b/pyproximal/proximal/L0.py index 4deecdf..f185f3b 100644 --- a/pyproximal/proximal/L0.py +++ b/pyproximal/proximal/L0.py @@ -3,6 +3,92 @@ from pyproximal.ProxOperator import _check_tau from pyproximal.projection import L0BallProj from pyproximal import ProxOperator +from pyproximal.proximal.L1 import _current_sigma + + +def _hardthreshold(x, thresh): + r"""Hard thresholding. + + Applies hard thresholding to vector ``x`` (equal to the proximity + operator for :math:`\|\mathbf{x}\|_0`) as shown in [1]_. + + .. [1] Chen, F., Shen, L., Suter, B.W., "Computing the proximity + operator of the Lp norm with 0 < p < 1", + IET Signal Processing, 10, 2016. + + Parameters + ---------- + x : :obj:`numpy.ndarray` + Vector + thresh : :obj:`float` + Threshold + + Returns + ------- + x1 : :obj:`numpy.ndarray` + Tresholded vector + + """ + x1 = x.copy() + x1[np.abs(x) <= thresh] = 0 + return x1 + + +class L0(ProxOperator): + r"""L0 norm proximal operator. + + Proximal operator of the :math:`\ell_0` norm: + :math:`\sigma\|\mathbf{x}\|_0 = \text{count}(x_i \ne 0)`. + + Parameters + ---------- + 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. + + Notes + ----- + The :math:`\ell_0` proximal operator is defined as: + + .. math:: + + \prox_{\tau \sigma \|\cdot\|_0}(\mathbf{x}) = + \operatorname{hard}(\mathbf{x}, \tau \sigma) = + \begin{cases} + x_i, & x_i < -\tau \sigma \\ + 0, & -\tau\sigma \leq x_i \leq \tau\sigma \\ + x_i, & x_i > \tau\sigma\\ + \end{cases} + + where :math:`\operatorname{hard}` is the so-called called *hard thresholding*. + + """ + def __init__(self, sigma=1.): + super().__init__(None, False) + self.sigma = sigma + self.count = 0 + + def __call__(self, x): + sigma = _current_sigma(self.sigma, self.count) + return np.sum(np.abs(x) > sigma) + + 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) + x = _hardthreshold(x, tau * sigma) + return x class L0Ball(ProxOperator): diff --git a/pyproximal/proximal/L1.py b/pyproximal/proximal/L1.py index ca0f84b..d1efb1e 100644 --- a/pyproximal/proximal/L1.py +++ b/pyproximal/proximal/L1.py @@ -67,7 +67,7 @@ class L1(ProxOperator): \operatorname{soft}(\mathbf{x}, \tau \sigma) = \begin{cases} x_i + \tau \sigma, & x_i - g_i < -\tau \sigma \\ - g_i, & -\sigma \leq x_i - g_i \leq \tau\sigma \\ + g_i, & -\tau\sigma \leq x_i - g_i \leq \tau\sigma \\ x_i - \tau\sigma, & x_i - g_i > \tau\sigma\\ \end{cases} diff --git a/pyproximal/proximal/__init__.py b/pyproximal/proximal/__init__.py index 5fed87b..fa933a4 100644 --- a/pyproximal/proximal/__init__.py +++ b/pyproximal/proximal/__init__.py @@ -10,6 +10,7 @@ AffineSet Affines set indicator Quadratic Quadratic function Nonlinear Nonlinear function + L0 L0 Norm L0Ball L0 Ball L1 L1 Norm L1Ball L1 Ball @@ -57,7 +58,7 @@ from .SingularValuePenalty import * __all__ = ['Box', 'Simplex', 'Intersection', 'AffineSet', 'Quadratic', - 'Euclidean', 'EuclideanBall', 'L0Ball', 'L1', 'L1Ball', 'L2', + 'Euclidean', 'EuclideanBall', 'L0', 'L0Ball', 'L1', 'L1Ball', 'L2', 'L2Convolve', 'L21', 'L21_plus_L1', 'Huber', 'Nuclear', 'NuclearBall', 'Orthogonal', 'VStack', 'Nonlinear', 'SCAD', 'Log', 'ETP', 'Geman', 'QuadraticEnvelopeCard', 'SingularValuePenalty']