# Source code for pyproximal.optimization.palm

import time
import numpy as np

def _backtracking(x, tau, H, proxf, ix, beta=0.5, niterback=10):
r"""Backtracking

Line-search algorithm for finding step sizes in palm algorithms when
the Lipschitz constant of the operator is unknown (or expensive to
estimate).

"""
def ftilde(x, y, f, g, tau, ix):
xy = x - y[ix]
return f(*y) + np.dot(g, xy) + \
(1. / (2. * tau)) * np.linalg.norm(xy) ** 2

iiterback = 0
if ix == 0:
else:
z = [x_.copy() for x_ in x]
while iiterback < niterback:
z[ix] = x[ix] - tau * grad
if proxf is not None:
z[ix] = proxf.prox(z[ix], tau)
ft = ftilde(z[ix], x, H, grad, tau, ix)
f = H(*z)
if f <= ft or tau < 1e-12:
break
tau *= beta
iiterback += 1
return z[ix], tau

[docs]def PALM(H, proxf, proxg, x0, y0, gammaf=1., gammag=1., beta=0.5, niter=10, niterback=100, callback=None, show=False): r"""Proximal Alternating Linearized Minimization Solves the following minimization problem using the Proximal Alternating Linearized Minimization (PALM) algorithm: .. math:: \mathbf{x}\mathbf{,y} = \argmin_{\mathbf{x}, \mathbf{y}} f(\mathbf{x}) + g(\mathbf{y}) + H(\mathbf{x}, \mathbf{y}) where :math:f(\mathbf{x}) and :math:g(\mathbf{y}) are any pair of convex functions that have known proximal operators, and :math:H(\mathbf{x}, \mathbf{y}) is a smooth function. Parameters ---------- H : :obj:pyproximal.utils.bilinear.Bilinear Bilinear function proxf : :obj:pyproximal.ProxOperator Proximal operator of f function proxg : :obj:pyproximal.ProxOperator Proximal operator of g function x0 : :obj:numpy.ndarray Initial x vector y0 : :obj:numpy.ndarray Initial y vector gammaf : :obj:float, optional Positive scalar weight for f function update. If None, use backtracking gammag : :obj:float, optional Positive scalar weight for g function update. If None, use backtracking beta : :obj:float, optional Backtracking parameter (must be between 0 and 1) niter : :obj:int, optional Number of iterations of iterative scheme niterback : :obj:int, optional Max number of iterations of backtracking callback : :obj:callable, optional Function with signature (callback(x)) to call after each iteration where x and y are the current model vectors show : :obj:bool, optional Display iterations log Returns ------- x : :obj:numpy.ndarray Inverted x vector y : :obj:numpy.ndarray Inverted y vector Notes ----- PALM [1]_ can be expressed by the following recursion: .. math:: \mathbf{x}^{k+1} = \prox_{c_k f}(\mathbf{x}^{k} - \frac{1}{c_k}\nabla_x H(\mathbf{x}^{k}, \mathbf{y}^{k}))\\ \mathbf{y}^{k+1} = \prox_{d_k g}(\mathbf{y}^{k} - \frac{1}{d_k}\nabla_y H(\mathbf{x}^{k+1}, \mathbf{y}^{k}))\\ Here :math:c_k=\gamma_f L_x and :math:d_k=\gamma_g L_y, where :math:L_x and :math:L_y are the Lipschitz constant of :math:\nabla_x H and :math:\nabla_y H, respectively. When such constants cannot be easily computed, a back-tracking algorithm can be instead employed to find suitable :math:c_k and :math:d_k parameters. .. [1] Bolte, J., Sabach, S., and Teboulle, M. "Proximal alternating linearized minimization for nonconvex and nonsmooth problems", Mathematical Programming, vol. 146, pp. 459–494. 2014. """ if show: tstart = time.time() print('PALM algorithm\n' '---------------------------------------------------------\n' 'Bilinear operator: %s\n' 'Proximal operator (f): %s\n' 'Proximal operator (g): %s\n' 'gammaf = %s\tgammag = %s\tniter = %d\n' % (type(H), type(proxf), type(proxg), str(gammaf), str(gammag), niter)) head = ' Itn x[0] y[0] f g H ck dk' print(head) backtrackingf, backtrackingg = False, False if gammaf is None: backtrackingf = True tauf = 1. ck = 0. if gammaf is None: backtrackingg = True taug = 1. dk = 0. x, y = x0.copy(), y0.copy() for iiter in range(niter): # x step if not backtrackingf: ck = gammaf * H.ly(y) x = x - (1. / ck) * H.gradx(x) if proxf is not None: x = proxf.prox(x, 1. / ck) else: x, tauf = _backtracking([x, y], tauf, H, proxf, 0, beta=beta, niterback=niterback) # update x parameter in H function H.updatex(x.copy()) # y step if not backtrackingg: dk = gammag * H.lx(x) y = y - (1. / dk) * H.grady(y) if proxg is not None: y = proxg.prox(y, 1. / dk) else: y, taug = _backtracking([x, y], tauf, H, proxf, 1, beta=beta, niterback=niterback) # update y parameter in H function 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 %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)) print('---------------------------------------------------------\n') return x, y
[docs]def iPALM(H, proxf, proxg, x0, y0, gammaf=1., gammag=1., a=[1., 1.], b=None, beta=0.5, niter=10, niterback=100, callback=None, show=False): r"""Inertial Proximal Alternating Linearized Minimization Solves the following minimization problem using the Inertial Proximal Alternating Linearized Minimization (iPALM) algorithm: .. math:: \mathbf{x}\mathbf{,y} = \argmin_{\mathbf{x}, \mathbf{y}} f(\mathbf{x}) + g(\mathbf{y}) + H(\mathbf{x}, \mathbf{y}) where :math:f(\mathbf{x}) and :math:g(\mathbf{y}) are any pair of convex functions that have known proximal operators, and :math:H(\mathbf{x}, \mathbf{y}) is a smooth function. Parameters ---------- H : :obj:pyproximal.utils.bilinear.Bilinear Bilinear function proxf : :obj:pyproximal.ProxOperator Proximal operator of f function proxg : :obj:pyproximal.ProxOperator Proximal operator of g function x0 : :obj:numpy.ndarray Initial x vector y0 : :obj:numpy.ndarray Initial y vector gammaf : :obj:float, optional Positive scalar weight for f function update. If None, use backtracking gammag : :obj:float, optional Positive scalar weight for g function update. If None, use backtracking a : :obj:list, optional Inertial parameters (:math:a \in [0, 1]) beta : :obj:float, optional Backtracking parameter (must be between 0 and 1) niter : :obj:int, optional Number of iterations of iterative scheme niterback : :obj:int, optional Max number of iterations of backtracking callback : :obj:callable, optional Function with signature (callback(x)) to call after each iteration where x and y are the current model vectors show : :obj:bool, optional Display iterations log Returns ------- x : :obj:numpy.ndarray Inverted x vector y : :obj:numpy.ndarray Inverted y vector Notes ----- iPALM [1]_ can be expressed by the following recursion: .. math:: \mathbf{x}_z^k = \mathbf{x}^k + \alpha_x (\mathbf{x}^k - \mathbf{x}^{k-1})\\ \mathbf{x}^{k+1} = \prox_{c_k f}(\mathbf{x}_z^k - \frac{1}{c_k}\nabla_x H(\mathbf{x}_z^k, \mathbf{y}^{k}))\\ \mathbf{y}_z^k = \mathbf{y}^k + \alpha_y (\mathbf{y}^k - \mathbf{y}^{k-1})\\ \mathbf{y}^{k+1} = \prox_{d_k g}(\mathbf{y}_z^k - \frac{1}{d_k}\nabla_y H(\mathbf{x}^{k+1}, \mathbf{y}_z^k)) Here :math:c_k=\gamma_f L_x and :math:d_k=\gamma_g L_y, where :math:L_x and :math:L_y are the Lipschitz constant of :math:\nabla_x H and :math:\nabla_y H, respectively. When such constants cannot be easily computed, a back-tracking algorithm can be instead employed to find suitable :math:c_k and :math:d_k parameters. Finally, note that we have implemented the version of iPALM where :math:\beta_x=\alpha_x and :math:\beta_y=\alpha_y. .. [1] Pock, T., and Sabach, S. "Inertial Proximal Alternating Linearized Minimization (iPALM) for Nonconvex and Nonsmooth Problems", SIAM Journal on Imaging Sciences, vol. 9. 2016. """ if show: tstart = time.time() print('iPALM algorithm\n' '---------------------------------------------------------\n' 'Bilinear operator: %s\n' 'Proximal operator (f): %s\n' 'Proximal operator (g): %s\n' 'gammaf = %s\tgammag = %s\n' 'a = %s\tniter = %d\n' % (type(H), type(proxf), type(proxg), str(gammaf), str(gammag), str(a), niter)) head = ' Itn x[0] y[0] f g H ck dk' print(head) backtrackingf, backtrackingg = False, False if gammaf is None: backtrackingf = True tauf = 1. ck = 0. if gammaf is None: backtrackingg = True taug = 1. dk = 0. x, y = x0.copy(), y0.copy() xold, yold = x0.copy(), y0.copy() for iiter in range(niter): # x step z = x + a[0] * (x - xold) if not backtrackingf: ck = gammaf * H.ly(y) xold = x.copy() x = z - (1. / ck) * H.gradx(z) if proxf is not None: x = proxf.prox(x, 1. / ck) else: xold = x.copy() x, tauf = _backtracking([z, y], tauf, H, proxf, 0, beta=beta, niterback=niterback) # update x parameter in H function H.updatex(x.copy()) # y step z = y + a[1] * (y - yold) if not backtrackingg: dk = gammag * H.lx(x) yold = y.copy() y = z - (1. / dk) * H.grady(z) if proxg is not None: y = proxg.prox(y, 1. / dk) else: yold = y.copy() y, taug = _backtracking([x, z], tauf, H, proxf, 1, beta=beta, niterback=niterback) # update y parameter in H function 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 %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)) print('---------------------------------------------------------\n') return x, y