PyUNLocBoX: Optimization by Proximal Splitting

The PyUNLocBoX is a Python package which uses proximal splitting methods to solve non-differentiable convex optimization problems. The documentation is available on Read the Docs and development takes place on GitHub. A (mostly unmaintained) Matlab version exists.

The package is designed to be easy to use while allowing any advanced tasks. It is not meant to be a black-box optimization tool. You’ll have to carefully design your solver. In exchange you’ll get full control of what the package does for you, without the pain of rewriting the proximity operators and the solvers and with the added benefit of tested algorithms. With this package, you can focus on your problem and the best way to solve it rather that the details of the algorithms.

Content

The following solvers are included:

• Gradient descent

• Forward-backward proximal splitting (FISTA and ISTA)

• Generalized forward-backward proximal splitting

• Douglas-Rachford proximal splitting

• Monotone+Lipschitz forward-backward-forward primal-dual

• Projection-based primal-dual

The following acceleration schemes are included:

• Backtracking acceleration based on a quadratic approximation of the objective

• FISTA acceleration for forward-backward solvers

• FISTA acceleration with backtracking for forward-backward solvers

• Regularized nonlinear acceleration (RNA) for gradient descent

To compose your objective, the following functions are included:

• L1-norm (eval, prox)

• L2-norm (eval, prox, grad)

• Nuclear-norm (eval, prox)

• TV-norm (eval, prox)

• Projection on the positive octant (eval, prox)

• Projection on the L2-ball (eval, prox)

• Structured sparsity (eval, prox)

Alternatively, you can easily define a custom function by implementing an evaluation method and a proximal operator or gradient method:

>>> from pyunlocbox import functions
>>> class myfunc(functions.func):
...     def _eval(self, x):
...         return 0  # Function evaluated at x.
...     def _grad(self, x):
...         return x  # Gradient evaluated at x, if available.
...     def _prox(self, x, T):
...         return x  # Proximal operator evaluated at x, if available.

Likewise, custom solvers are defined by inheriting from solvers.solver and implementing _pre, _algo, and _post. Custom acceleration schemes are defined by inheriting from acceleration.accel and implementing _pre, _update_step, _update_sol, and _post.

Usage

Following is a typical usage example that solves an optimization problem composed by the sum of two convex functions. The functions and solver objects are first instantiated with the desired parameters. The problem is then solved by a call to the solving function.

>>> from pyunlocbox import functions, solvers
>>> f1 = functions.norm_l2(y=[4, 5, 6, 7])
>>> f2 = functions.dummy()
>>> solver = solvers.forward_backward()
>>> ret = solvers.solve([f1, f2], [0., 0, 0, 0], solver, atol=1e-5)
Solution found after 9 iterations:
    objective function f(sol) = 6.714385e-08
    stopping criterion: ATOL
>>> ret['sol']
array([3.99990766, 4.99988458, 5.99986149, 6.99983841])

You can try it online, look at the tutorials to learn how to use it, or look at the reference guide for an exhaustive documentation of the API. Enjoy!

Installation

The PyUNLocBoX is available on PyPI:

$ pip install pyunlocbox

The PyUNLocBoX is available on conda-forge:

$ conda install -c conda-forge pyunlocbox

Contributing

See the guidelines for contributing in CONTRIBUTING.rst.

Similar libraries

Other proximal based algorithms and operators can be found in:

• http://proximity-operator.net/

• https://gitlab.gwdg.de/nam/ProxPython

Furthermore, many proximal operators are availlable in the proxop python library.

Acknowledgments

The PyUNLocBoX was started in 2014 as an academic open-source project for research purpose at the EPFL LTS2 laboratory.

It is released under the terms of the BSD 3-Clause license.

If you are using the library for your research, for the sake of reproducibility, please cite the version you used as indexed by Zenodo. Or cite the generic concept as:

@misc{pyunlocbox,
  title = {PyUNLocBoX: Optimization by Proximal Splitting},
  author = {Defferrard, Micha\"el and Pena, Rodrigo and Perraudin, Nathana\"el},
  doi = {10.5281/zenodo.1199081},
  url = {https://github.com/epfl-lts2/pyunlocbox/},
}

Tutorials

The following are some tutorials which show and explain how to use the toolbox to solve some real problems. They goes in increasing degree of difficulty. If you have never used the toolbox before, you are encouraged to follow them in order as they build one upon the other.

Simple least square problem

This simplistic example is only meant to demonstrate the basic workflow of the toolbox. Here we want to solve a least square problem, i.e. we want the solution to converge to the original signal without any constraint. Lets define this signal by :

>>> y = [4, 5, 6, 7]

The first function to minimize is the sum of squared distances between the current signal x and the original y. For this purpose, we instantiate an L2-norm object :

>>> from pyunlocbox import functions
>>> f1 = functions.norm_l2(y=y)

This standard function object provides the eval(), grad() and prox() methods that will be useful to the solver. We can evaluate them at any given point :

>>> f1.eval([0, 0, 0, 0])
126
>>> f1.grad([0, 0, 0, 0])
array([ -8, -10, -12, -14])
>>> f1.prox([0, 0, 0, 0], 1)
array([2.66666667, 3.33333333, 4.        , 4.66666667])

We need a second function to minimize, which usually describes a constraint. As we have no constraint, we just define a dummy function object by hand. We have to define the _eval() and _grad() methods as the solver we will use requires it :

>>> f2 = functions.func()
>>> f2._eval = lambda x: 0
>>> f2._grad = lambda x: 0

Note

We could also have used the pyunlocbox.functions.dummy function object.

We can now instantiate the solver object :

>>> from pyunlocbox import solvers
>>> solver = solvers.forward_backward()

And finally solve the problem :

>>> x0 = [0., 0., 0., 0.]
>>> ret = solvers.solve([f2, f1], x0, solver, atol=1e-5, verbosity='HIGH')
INFO: Forward-backward method
    func evaluation: 0.000000e+00
    norm_l2 evaluation: 1.260000e+02
Iteration 1 of forward_backward:
    func evaluation: 0.000000e+00
    norm_l2 evaluation: 1.400000e+01
    objective = 1.40e+01
Iteration 2 of forward_backward:
    func evaluation: 0.000000e+00
    norm_l2 evaluation: 2.963739e-01
    objective = 2.96e-01
Iteration 3 of forward_backward:
    func evaluation: 0.000000e+00
    norm_l2 evaluation: 7.902529e-02
    objective = 7.90e-02
Iteration 4 of forward_backward:
    func evaluation: 0.000000e+00
    norm_l2 evaluation: 5.752265e-02
    objective = 5.75e-02
Iteration 5 of forward_backward:
    func evaluation: 0.000000e+00
    norm_l2 evaluation: 5.142032e-03
    objective = 5.14e-03
Iteration 6 of forward_backward:
    func evaluation: 0.000000e+00
    norm_l2 evaluation: 1.553851e-04
    objective = 1.55e-04
Iteration 7 of forward_backward:
    func evaluation: 0.000000e+00
    norm_l2 evaluation: 5.498523e-04
    objective = 5.50e-04
Iteration 8 of forward_backward:
    func evaluation: 0.000000e+00
    norm_l2 evaluation: 1.091372e-04
    objective = 1.09e-04
Iteration 9 of forward_backward:
    func evaluation: 0.000000e+00
    norm_l2 evaluation: 6.714385e-08
    objective = 6.71e-08
Solution found after 9 iterations:
    objective function f(sol) = 6.714385e-08
    stopping criterion: ATOL

The solving function returns several values, one is the found solution :

>>> ret['sol']
array([3.99990766, 4.99988458, 5.99986149, 6.99983841])

Another one is the value returned by each function objects at each iteration. As we passed two function objects (L2-norm and dummy), the objective is a 2 by 11 (10 iterations plus the evaluation at x0) ndarray. Lets plot a convergence graph out of it :

>>> import numpy as np
>>> objective = np.array(ret['objective'])
>>> import matplotlib.pyplot as plt
>>> _ = plt.figure()
>>> _ = plt.semilogy(objective[:, 1], 'x', label='L2-norm')
>>> _ = plt.semilogy(objective[:, 0], label='Dummy')
>>> _ = plt.semilogy(np.sum(objective, axis=1), label='Global objective')
>>> _ = plt.grid(True)
>>> _ = plt.title('Convergence')
>>> _ = plt.legend(numpoints=1)
>>> _ = plt.xlabel('Iteration number')
>>> _ = plt.ylabel('Objective function value')

The above graph shows an exponential convergence of the objective function. The global objective is obviously only composed of the L2-norm as the dummy function object was defined to always evaluate to 0 (f2._eval = lambda x: 0).

Compressed sensing using forward-backward

This tutorial presents a compressed sensing problem solved by the forward-backward splitting algorithm. The convex optimization problem is the sum of a data fidelity term and a regularization term which expresses a prior on the sparsity of the solution, given by

\[\min\limits_x \|Ax-y\|_2^2 + \tau \|x\|_1\]

where y are the measurements, A is the measurement matrix and \(\tau\) expresses the trade-off between the two terms.

The number of necessary measurements m is computed with respect to the signal size n and the sparsity level S in order to very often perform a perfect reconstruction. See [CR07] for details.

>>> n = 5000
>>> S = 100
>>> import numpy as np
>>> m = int(np.ceil(S * np.log(n)))
>>> print('Number of measurements: {}'.format(m))
Number of measurements: 852
>>> print('Compression ratio: {:3.2f}'.format(float(n) / m))
Compression ratio: 5.87

We generate a random measurement matrix A:

>>> np.random.seed(1)  # Reproducible results.
>>> A = np.random.normal(size=(m, n))

Create the S sparse signal x:

>>> x = np.zeros(n)
>>> I = np.random.permutation(n)
>>> x[I[0:S]] = np.random.normal(size=S)
>>> x = x / np.linalg.norm(x)

Generate the measured signal y:

>>> y = np.dot(A, x)

The prior objective to minimize is defined by

\[f_1(x) = \tau \|x\|_1\]

which can be expressed by the toolbox L1-norm function object. It can be instantiated as follows, while setting the regularization parameter tau:

>>> from pyunlocbox import functions
>>> tau = 1.0
>>> f1 = functions.norm_l1(lambda_=tau)

The fidelity objective to minimize is defined by

\[f_2(x) = \|Ax-y\|_2^2\]

which can be expressed by the toolbox L2-norm function object. It can be instantiated as follows:

>>> f2 = functions.norm_l2(y=y, A=A)

or alternatively as follows:

>>> f3 = functions.norm_l2(y=y)
>>> f3.A = lambda x: np.dot(A, x)
>>> f3.At = lambda x: np.dot(A.T, x)

Note

In this case the forward and adjoint operators were passed as functions not as matrices.

A third alternative would be to define the function object by hand:

>>> f4 = functions.func()
>>> f4._grad = lambda x: 2.0 * np.dot(A.T, np.dot(A, x) - y)
>>> f4._eval = lambda x: np.linalg.norm(np.dot(A, x) - y)**2

Note

The three alternatives to instantiate the function objects (f2, f3 and f4) are strictly equivalent and give the exact same results.

Now that the two function objects to minimize (the L1-norm and the L2-norm) are instantiated, we can instantiate the solver object. The step size for optimal convergence is \(\frac{1}{\beta}\) where \(\beta\) is the Lipschitz constant of the gradient of f2, f3, f4 given by:

\[\beta = 2 \cdot \|A\|_{\text{op}}^2 = 2 \cdot \lambda_{max} (A^*A).\]

To solve this problem, we use the Forward-Backward splitting algorithm which is instantiated as follows:

>>> step = 0.5 / np.linalg.norm(A, ord=2)**2
>>> from pyunlocbox import solvers
>>> solver = solvers.forward_backward(step=step)

Note

A complete description of the constructor parameters and default values is given by the solver object pyunlocbox.solvers.forward_backward reference documentation.

After the instantiations of the functions and solver objects, the setting of a starting point x0, the problem is solved by the toolbox solving function as follows:

>>> x0 = np.zeros(n)
>>> ret = solvers.solve([f1, f2], x0, solver, rtol=1e-4, maxit=300)
Solution found after 151 iterations:
    objective function f(sol) = 7.668167e+00
    stopping criterion: RTOL

Note

A complete description of the parameters, their default values and the returned values is given by the solving function pyunlocbox.solvers.solve() reference documentation.

Let’s display the results:

>>> import matplotlib.pyplot as plt
>>> _ = plt.figure()
>>> _ = plt.plot(x, 'o', label='Original')
>>> _ = plt.plot(ret['sol'], 'xr', label='Reconstructed')
>>> _ = plt.grid(True)
>>> _ = plt.title('Achieved reconstruction')
>>> _ = plt.legend(numpoints=1)
>>> _ = plt.xlabel('Signal dimension number')
>>> _ = plt.ylabel('Signal value')

The above figure shows a good reconstruction which is both sparse (thanks to the L1-norm objective) and close to the measurements (thanks to the L2-norm objective).

Let’s display the convergence of the two objective functions:

>>> objective = np.array(ret['objective'])
>>> _ = plt.figure()
>>> _ = plt.semilogy(objective[:, 0], label='L1-norm objective')
>>> _ = plt.semilogy(objective[:, 1], label='L2-norm objective')
>>> _ = plt.semilogy(np.sum(objective, axis=1), label='Global objective')
>>> _ = plt.grid(True)
>>> _ = plt.title('Convergence')
>>> _ = plt.legend()
>>> _ = plt.xlabel('Iteration number')
>>> _ = plt.ylabel('Objective function value')

Compressed sensing using Douglas-Rachford

This tutorial presents a compressed sensing problem solved by the Douglas-Rachford splitting algorithm. The convex optimization problem, a term which expresses a prior on the sparsity of the solution constrained by some data fidelity, is given by

\[\min\limits_x \|x\|_1 \text{ s.t. } \|Ax-y\|_2 \leq \epsilon\]

where y are the measurements and A is the measurement matrix.

The number of necessary measurements m is computed with respect to the signal size n and the sparsity level S in order to very often perform a perfect reconstruction. See [CR07] for details.

>>> n = 900
>>> S = 45
>>> import numpy as np
>>> m = int(np.ceil(S * np.log(n)))
>>> print('Number of measurements: {}'.format(m))
Number of measurements: 307
>>> print('Compression ratio: {:3.2f}'.format(float(n) / m))
Compression ratio: 2.93

We generate a random measurement matrix A:

>>> np.random.seed(1)  # Reproducible results.
>>> A = np.random.normal(size=(m, n))

Create the S sparse signal x:

>>> x = np.zeros(n)
>>> I = np.random.permutation(n)
>>> x[I[0:S]] = np.random.normal(size=S)
>>> x = x / np.linalg.norm(x)

Generate the measured signal y:

>>> y = np.dot(A, x)

The first objective function to minimize is defined by

\[f_1(x) = \|x\|_1\]

which can be expressed by the toolbox L1-norm function object. It can be instantiated as follows:

>>> from pyunlocbox import functions
>>> f1 = functions.norm_l1()

The second objective function to minimize is defined by

\[f_2(x) = \iota_C(x)\]

where \(\iota_C()\) is the indicator function of the set \(C = \left\{z \in \mathbb{R}^n \mid \|Az-y\|_2 \leq \epsilon \right\}\) which is zero if \(z\) is in the set and infinite otherwise. This function can be expressed by the toolbox L2-ball function object which can be instantiated as follows:

>>> f2 = functions.proj_b2(epsilon=1e-7, y=y, A=A, tight=False,
... nu=np.linalg.norm(A, ord=2)**2)

Now that the two function objects to minimize (the L1-norm and the L2-ball) are instantiated, we can instantiate the solver object. To solve this problem, we use the Douglas-Rachford splitting algorithm which is instantiated as follows:

>>> from pyunlocbox import solvers
>>> solver = solvers.douglas_rachford(step=1e-2)

After the instantiations of the functions and solver objects, the setting of a starting point x0, the problem is solved by the toolbox solving function as follows:

>>> x0 = np.zeros(n)
>>> ret = solvers.solve([f1, f2], x0, solver, rtol=1e-4, maxit=300)
Solution found after 43 iterations:
    objective function f(sol) = 5.607407e+00
    stopping criterion: RTOL

Let’s display the results:

>>> import matplotlib.pyplot as plt
>>> _ = plt.figure()
>>> _ = plt.plot(x, 'o', label='Original')
>>> _ = plt.plot(ret['sol'], 'xr', label='Reconstructed')
>>> _ = plt.grid(True)
>>> _ = plt.title('Achieved reconstruction')
>>> _ = plt.legend(numpoints=1)
>>> _ = plt.xlabel('Signal dimension number')
>>> _ = plt.ylabel('Signal value')

The above figure shows a good reconstruction which is both sparse (thanks to the L1-norm objective) and close to the measurements (thanks to the L2-ball constraint).

Let’s display the convergence of the objective function:

>>> objective = np.array(ret['objective'])
>>> _ = plt.figure()
>>> _ = plt.semilogy(objective[:, 0], label='L1-norm objective')
>>> _ = plt.grid(True)
>>> _ = plt.title('Convergence')
>>> _ = plt.legend()
>>> _ = plt.xlabel('Iteration number')
>>> _ = plt.ylabel('Objective function value')

Image reconstruction (Forward-Backward, Total Variation, L2-norm)

This tutorial presents an image reconstruction problem solved by the Forward-Backward splitting algorithm. The convex optimization problem is the sum of a data fidelity term and a regularization term which expresses a prior on the smoothness of the solution, given by

\[\min\limits_x \tau \|g(x)-y\|_2^2 + \|x\|_\text{TV}\]

where \(\|\cdot\|_\text{TV}\) denotes the total variation, y are the measurements, g is a masking operator and \(\tau\) expresses the trade-off between the two terms.

Load an image and convert it to grayscale

>>> import matplotlib.image as mpimg
>>> import numpy as np
>>> try:
...     im_original = mpimg.imread('tutorials/barbara.png')
... except:
...     im_original = mpimg.imread('doc/tutorials/barbara.png')

and generate a random masking matrix

>>> np.random.seed(14)  # Reproducible results.
>>> mask = np.random.uniform(size=im_original.shape)
>>> mask = mask > 0.85

which masks 85% of the pixels. The masked image is given by

>>> g = lambda x: mask * x
>>> im_masked = g(im_original)

The prior objective to minimize is defined by

\[f_1(x) = \|x\|_\text{TV}\]

which can be expressed by the toolbox TV-norm function object, instantiated with

>>> from pyunlocbox import functions
>>> f1 = functions.norm_tv(maxit=50, dim=2)

The fidelity objective to minimize is defined by

\[f_2(x) = \tau \|g(x)-y\|_2^2\]

which can be expressed by the toolbox L2-norm function object, instantiated with

>>> tau = 100
>>> f2 = functions.norm_l2(y=im_masked, A=g, lambda_=tau)

Note

We set \(\tau\) to a large value as we trust our measurements and want the solution to be close to them. For noisy measurements a lower value should be considered.

The step size for optimal convergence is \(\frac{1}{\beta}\) where \(\beta=2\tau\) is the Lipschitz constant of the gradient of \(f_2\) [BT09a]. The Forward-Backward splitting algorithm is instantiated with

>>> from pyunlocbox import solvers
>>> solver = solvers.forward_backward(step=0.5/tau)

and the problem solved with

>>> x0 = np.array(im_masked)  # Make a copy to preserve im_masked.
>>> ret = solvers.solve([f1, f2], x0, solver, maxit=100)
Solution found after 78 iterations:
    objective function f(sol) = 6.723857e+03
    stopping criterion: RTOL

Let’s display the results:

>>> import matplotlib.pyplot as plt
>>> fig = plt.figure(figsize=(8, 2.5))
>>> ax1 = fig.add_subplot(1, 3, 1)
>>> _ = ax1.imshow(im_original, cmap='gray')
>>> _ = ax1.axis('off')
>>> _ = ax1.set_title('Original image')
>>> ax2 = fig.add_subplot(1, 3, 2)
>>> _ = ax2.imshow(im_masked, cmap='gray')
>>> _ = ax2.axis('off')
>>> _ = ax2.set_title('Masked image')
>>> ax3 = fig.add_subplot(1, 3, 3)
>>> _ = ax3.imshow(ret['sol'], cmap='gray')
>>> _ = ax3.axis('off')
>>> _ = ax3.set_title('Reconstructed image')

The above figure shows a good reconstruction which is both smooth (the TV prior) and close to the measurements (the L2 fidelity).

Image denoising (Douglas-Rachford, Total Variation, L2-norm)

This tutorial presents an image denoising problem solved by the Douglas-Rachford splitting algorithm. The convex optimization problem, a term which expresses a prior on the smoothness of the solution constrained by some data fidelity, is given by

\[\min\limits_x \|x\|_\text{TV} \text{ s.t. } \|x-y\|_2 \leq \epsilon\]

where \(\|\cdot\|_\text{TV}\) denotes the total variation, y are the measurements and \(\epsilon\) expresses the noise level.

Create a white circle on a black background

>>> import numpy as np
>>> N = 650
>>> im_original = np.resize(np.linspace(-1, 1, N), (N, N))
>>> im_original = np.sqrt(im_original**2 + im_original.T**2)
>>> im_original = im_original < 0.7

and add some random Gaussian noise

>>> sigma = 0.5  # Variance of 0.25.
>>> np.random.seed(7)  # Reproducible results.
>>> im_noisy = im_original + sigma * np.random.normal(size=im_original.shape)

The prior objective function to minimize is defined by

\[f_1(x) = \|x\|_\text{TV}\]

which can be expressed by the toolbox TV-norm function object, instantiated with

>>> from pyunlocbox import functions
>>> f1 = functions.norm_tv(maxit=50, dim=2)

The fidelity constraint expressed as an objective function to minimize is defined by

\[f_2(x) = \iota_S(x)\]

where \(\iota_S()\) is the indicator function of the set \(S = \left\{z \in \mathbb{R}^n \mid \|z-y\|_2 \leq \epsilon \right\}\) which is zero if \(z\) is in the set and infinite otherwise. This function can be expressed by the toolbox L2-ball function, instantiated with

>>> y = np.reshape(im_noisy, -1)  # Reshape the 2D image as a 1D vector.
>>> epsilon = N * sigma           # Variance multiplied by N^2.
>>> f = functions.proj_b2(y=y, epsilon=epsilon)
>>> f2 = functions.func()
>>> f2._eval = lambda x: 0        # Indicator functions evaluate to zero.
>>> def prox(x, step):
...     return np.reshape(f.prox(np.reshape(x, -1), 0), im_noisy.shape)
>>> f2._prox = prox

Note

We defined a custom proximal operator which transforms the 2D image as a 1D vector because pyunlocbox.functions.proj_b2 operates on the columns of x while pyunlocbox.functions.norm_tv needs a two-dimensional array to compute the 2D TV norm.

The Douglas-Rachford splitting algorithm is instantiated with

>>> from pyunlocbox import solvers
>>> solver = solvers.douglas_rachford(step=0.1)

and the problem solved with

>>> x0 = np.array(im_noisy)  # Make a copy to preserve y aka im_noisy.
>>> ret = solvers.solve([f1, f2], x0, solver)
Solution found after 25 iterations:
    objective function f(sol) = 2.080376e+03
    stopping criterion: RTOL

Let’s display the results:

>>> import matplotlib.pyplot as plt
>>> fig = plt.figure(figsize=(8, 2.5))
>>> ax1 = fig.add_subplot(1, 3, 1)
>>> _ = ax1.imshow(im_original, cmap='gray')
>>> _ = ax1.axis('off')
>>> _ = ax1.set_title('Original image')
>>> ax2 = fig.add_subplot(1, 3, 2)
>>> _ = ax2.imshow(im_noisy, cmap='gray')
>>> _ = ax2.axis('off')
>>> _ = ax2.set_title('Noisy image')
>>> ax3 = fig.add_subplot(1, 3, 3)
>>> _ = ax3.imshow(ret['sol'], cmap='gray')
>>> _ = ax3.axis('off')
>>> _ = ax3.set_title('Denoised image')

The above figure shows a good reconstruction which is both smooth (the TV prior) and close to the measurements (the L2 fidelity constraint).

API reference

The package is mainly organized around two class hierarchies: the functions and the solvers. Instantiated functions represent convex functions to optimize. Instantiated solvers represent solving algorithms. The pyunlocbox.solvers.solve() solving function takes as parameters a solver object and some function objects to actually solve the optimization problem. See this function’s documentation for a typical usage example.

The pyunlocbox package is divided into the following modules:

• functions: objective functions to define an optimization problem,

• solvers: the main solving function and common solvers,

• acceleration: general acceleration schemes for various solvers,

• operators: some operators.

pyunlocbox.test()

Run the test suite.

Functions

The pyunlocbox.functions module implements an interface for solvers to access the functions to be optimized as well as common objective functions.

Interface

The func base class defines a common interface to all functions:

func.cap(x)	Test the capabilities of the function object.
func.eval(x)	Function evaluation.
func.prox(x, T)	Function proximal operator.
func.grad(x)	Function gradient.

Functions

Then, derived classes implement various common objective functions.

Norm operators (based on norm)

norm_l1(**kwargs)	L1-norm (eval, prox).
norm_l2(**kwargs)	L2-norm (eval, prox, grad).
norm_nuclear(**kwargs)	Nuclear-norm (eval, prox).
norm_tv([dim, verbosity])	TV-norm (eval, prox).

Projection operators (based on proj)

proj_positive(**kwargs)	Projection on the positive octant (eval, prox).
proj_b2([epsilon, method])	Projection on the L2-ball (eval, prox).
proj_lineq([A, pinvA])	Projection on the plane satisfying the linear equality Az = y (eval, prox).
proj_spsd(**kwargs)	Projection on symmetric positive semi-definite matrices (eval, prox).

Miscellaneous

dummy(**kwargs)	Dummy function (eval, prox, grad).
structured_sparsity([lambda_, groups, weights])	Structured sparsity (eval, prox).

class pyunlocbox.functions.func(y=0, A=None, At=None, tight=True, nu=1, tol=0.001, maxit=200, **kwargs)

Bases: object

This class defines the function object interface.

It is intended to be a base class for standard functions which will implement the required methods. It can also be instantiated by user code and dynamically modified for rapid testing. The instanced objects are meant to be passed to the pyunlocbox.solvers.solve() solving function.

Parameters

yarray_like, optional

Measurements. Default is 0.

Afunction or ndarray, optional

The forward operator. Default is the identity, \(A(x)=x\). If A is an ndarray, it will be converted to the operator form.

Atfunction or ndarray, optional

The adjoint operator. If At is an ndarray, it will be converted to the operator form. If A is an ndarray, default is the transpose of A. If A is a function, default is A, \(At(x)=A(x)\).

tightbool, optional

True if A is a tight frame (semi-orthogonal linear transform), False otherwise. Default is True.

nufloat, optional

Bound on the norm of the operator A, i.e. \(\|A(x)\|^2 \leq \nu \|x\|^2\). Default is 1.

tolfloat, optional

The tolerance stopping criterion. The exact definition depends on the function object, please see the documentation of the considered function. Default is 1e-3.

maxitint, optional

The maximum number of iterations. Default is 200.

Examples

Let’s define a parabola as an example of the manual implementation of a function object :

>>> from pyunlocbox import functions
>>> f = functions.func()
>>> f._eval = lambda x: x**2
>>> f._grad = lambda x: 2*x
>>> x = [1, 2, 3, 4]
>>> f.eval(x)
array([ 1,  4,  9, 16])
>>> f.grad(x)
array([2, 4, 6, 8])
>>> f.cap(x)
['EVAL', 'GRAD']

eval(x)

Function evaluation.

Parameters

xarray_like

The evaluation point. If x is a matrix, the function gets evaluated for each column, as if it was a set of independent problems. Some functions, like the nuclear norm, are only defined on matrices.

Returns

zfloat

The objective function evaluated at x. If x is a matrix, the sum of the objectives is returned.

Notes

This method is required by the pyunlocbox.solvers.solve() solving function to evaluate the objective function. Each function class should therefore define it.

prox(x, T)

Function proximal operator.

Parameters

xarray_like

The evaluation point. If x is a matrix, the function gets evaluated for each column, as if it was a set of independent problems. Some functions, like the nuclear norm, are only defined on matrices.

Tfloat

The regularization parameter.

Returns

zndarray

The proximal operator evaluated for each column of x.

Notes

The proximal operator is defined by \(\operatorname{prox}_{\gamma f}(x) = \operatorname{arg\,min} \limits_z \frac{1}{2} \|x-z\|_2^2 + \gamma f(z)\)

This method is required by some solvers.

When the map A in the function construction is a tight frame (semi-orthogonal linear transformation), we can use property (x) of Table 10.1 in [CP11] to compute the proximal operator of the composition of A with the base function. Whenever this is not the case, we have to resort to some iterative procedure, which may be very inefficient.

grad(x)

Function gradient.

Parameters

xarray_like

The evaluation point. If x is a matrix, the function gets evaluated for each column, as if it was a set of independent problems. Some functions, like the nuclear norm, are only defined on matrices.

Returns

zndarray

The objective function gradient evaluated for each column of x.

Notes

This method is required by some solvers.

cap(x)

Test the capabilities of the function object.

Parameters

xarray_like

The evaluation point. Not really needed, but this function calls the methods of the object to test if they can properly execute without raising an exception. Therefore it needs some evaluation point with a consistent size.

Returns

caplist of string

A list of capabilities (‘EVAL’, ‘GRAD’, ‘PROX’).

class pyunlocbox.functions.dummy(**kwargs)

Bases: pyunlocbox.functions.func

Dummy function (eval, prox, grad).

This can be used as a second function object when there is only one function to minimize. It always evaluates as 0.

Examples

>>> from pyunlocbox import functions
>>> f = functions.dummy()
>>> x = [1, 2, 3, 4]
>>> f.eval(x)
0
>>> f.prox(x, 1)
array([1, 2, 3, 4])
>>> f.grad(x)
array([0, 0, 0, 0])

class pyunlocbox.functions.norm(lambda_=1, w=1, **kwargs)

Bases: pyunlocbox.functions.func

Base class which defines the attributes of the norm objects.

See generic attributes descriptions of the pyunlocbox.functions.func base class.

Parameters

lambda_float, optional

Regularization parameter \(\lambda\). Default is 1.

warray_like, optional

Weights for a weighted norm. Default is 1.

class pyunlocbox.functions.norm_l1(**kwargs)

Bases: pyunlocbox.functions.norm

L1-norm (eval, prox).

See generic attributes descriptions of the pyunlocbox.functions.norm base class. Note that the constructor takes keyword-only parameters.

Notes

• The L1-norm of the vector x is given by \(\lambda \|w \cdot (A(x)-y)\|_1\).

• The L1-norm proximal operator evaluated at x is given by \(\operatorname{arg\,min}\limits_z \frac{1}{2} \|x-z\|_2^2 + \gamma \|w \cdot (A(z)-y)\|_1\) where \(\gamma = \lambda \cdot T\). This is simply a soft thresholding.

Examples

>>> from pyunlocbox import functions
>>> f = functions.norm_l1()
>>> f.eval([1, 2, 3, 4])
10
>>> f.prox([1, 2, 3, 4], 1)
array([0., 1., 2., 3.])

class pyunlocbox.functions.norm_l2(**kwargs)

Bases: pyunlocbox.functions.norm

L2-norm (eval, prox, grad).

See generic attributes descriptions of the pyunlocbox.functions.norm base class. Note that the constructor takes keyword-only parameters.

Notes

• The squared L2-norm of the vector x is given by \(\lambda \|w \cdot (A(x)-y)\|_2^2\).

• The squared L2-norm proximal operator evaluated at x is given by \(\operatorname{arg\,min}\limits_z \frac{1}{2} \|x-z\|_2^2 + \gamma \|w \cdot (A(z)-y)\|_2^2\) where \(\gamma = \lambda \cdot T\).

• The squared L2-norm gradient evaluated at x is given by \(2 \lambda \cdot At(w \cdot (A(x)-y))\).

Examples

>>> from pyunlocbox import functions
>>> f = functions.norm_l2()
>>> x = [1, 2, 3, 4]
>>> f.eval(x)
30
>>> f.prox(x, 1)
array([0.33333333, 0.66666667, 1.        , 1.33333333])
>>> f.grad(x)
array([2, 4, 6, 8])

class pyunlocbox.functions.norm_nuclear(**kwargs)

Bases: pyunlocbox.functions.norm

Nuclear-norm (eval, prox).

See generic attributes descriptions of the pyunlocbox.functions.norm base class. Note that the constructor takes keyword-only parameters.

Notes

• The nuclear-norm of the matrix x is given by \(\lambda \| x \|_* = \lambda \operatorname{trace} (\sqrt{x^* x}) = \lambda \sum_{i=1}^N |e_i|\) where e_i are the eigenvalues of x.

• The nuclear-norm proximal operator evaluated at x is given by \(\operatorname{arg\,min}\limits_z \frac{1}{2} \|x-z\|_2^2 + \gamma \| x \|_*\) where \(\gamma = \lambda \cdot T\), which is a soft-thresholding of the eigenvalues.

Examples

>>> from pyunlocbox import functions
>>> f = functions.norm_nuclear()
>>> f.eval([[1, 2],[2, 3]])  
4.47213595...
>>> f.prox([[1, 2],[2, 3]], 1)
array([[0.89442719, 1.4472136 ],
       [1.4472136 , 2.34164079]])

class pyunlocbox.functions.norm_tv(dim=2, verbosity='LOW', **kwargs)

Bases: pyunlocbox.functions.norm

TV-norm (eval, prox).

See generic attributes descriptions of the pyunlocbox.functions.norm base class. Note that the constructor takes keyword-only parameters.

Notes

TODO

See [BT09b] for details about the algorithm.

Examples

>>> import numpy as np
>>> from pyunlocbox import functions
>>> f = functions.norm_tv()
>>> x = np.arange(0, 16)
>>> x = x.reshape(4, 4)
>>> f.eval(x)  
    norm_tv evaluation: 5.210795e+01
52.10795063...

class pyunlocbox.functions.proj(**kwargs)

Bases: pyunlocbox.functions.func

Base class which defines the attributes of the proj objects.

See generic attributes descriptions of the pyunlocbox.functions.func base class.

Notes

• All indicator functions (projections) evaluate to zero by definition.

class pyunlocbox.functions.proj_positive(**kwargs)

Bases: pyunlocbox.functions.proj

Projection on the positive octant (eval, prox).

This function is the indicator function \(i_S(z)\) of the set \(S = \left\{z \in \mathbb{R}^N \mid z \leq 0 \right\}\) that is zero if \(z\) is in the set and infinite otherwise.

See generic attributes descriptions of the pyunlocbox.functions.proj base class. Note that the constructor takes keyword-only parameters.

Notes

• The evaluation of this function is zero.

Examples

>>> from pyunlocbox import functions
>>> f = functions.proj_positive()
>>> x = [-2.5, 1.5]
>>> f.eval(x)
0
>>> f.prox(x, 0)
array([0. , 1.5])

class pyunlocbox.functions.proj_spsd(**kwargs)

Bases: pyunlocbox.functions.proj

Projection on symmetric positive semi-definite matrices (eval, prox).

This function is the indicator function \(i_S(M)\) of the set \(S = \left\{M \in \mathbb{R}^{N \times N} \mid M \succeq 0, M=M^T \right\}\) that is zero if \(M\) is in the set and infinite otherwise.

See generic attributes descriptions of the pyunlocbox.functions.proj base class. Note that the constructor takes keyword-only parameters.

Notes

• The evaluation of this function is zero.

Examples

>>> from pyunlocbox import functions
>>> f = functions.proj_spsd()
>>> A = np.array([[0, -1] , [-1, 1]])
>>> A = (A + A.T) / 2  # Symmetrize the matrix.
>>> np.linalg.eig(A)[0]
array([-0.61803399,  1.61803399])
>>> f.eval(A)
0
>>> Aproj = f.prox(A, 0)
>>> np.linalg.eig(Aproj)[0]
array([0.        , 1.61803399])

class pyunlocbox.functions.proj_b2(epsilon=1, method='FISTA', **kwargs)

Bases: pyunlocbox.functions.proj

Projection on the L2-ball (eval, prox).

This function is the indicator function \(i_S(z)\) of the set \(S= \left\{z \in \mathbb{R}^N \mid \|Az-y\|_2 \leq \epsilon \right\}\) that is zero if \(z\) is in the set and infinite otherwise.

See generic attributes descriptions of the pyunlocbox.functions.proj base class. Note that the constructor takes keyword-only parameters.

Parameters

epsilonfloat, optional

The radius of the ball. Default is 1.

method{‘FISTA’, ‘ISTA’}, optional

The method used to solve the problem. It can be ‘FISTA’ or ‘ISTA’. Default is ‘FISTA’.

See also

proj_lineq

use instead of epsilon=0

Notes

• The tol parameter is defined as the tolerance for the projection on the L2-ball. The algorithm stops if \(\frac{\epsilon}{1-tol} \leq \|y-A(z)\|_2 \leq \frac{\epsilon}{1+tol}\).

• The evaluation of this function is zero.

• The L2-ball proximal operator evaluated at x is given by \(\operatorname{arg\,min}\limits_z \frac{1}{2} \|x-z\|_2^2 + i_S(z)\) which has an identical solution as \(\operatorname{arg\,min}\limits_z \|x-z\|_2^2\) such that \(\|A(z)-y\|_2 \leq \epsilon\). It is thus a projection of the vector x onto an L2-ball of diameter epsilon.

Examples

>>> from pyunlocbox import functions
>>> f = functions.proj_b2(y=[1, 1])
>>> x = [3, 3]
>>> f.eval(x)
0
>>> f.prox(x, 0)
array([1.70710678, 1.70710678])

class pyunlocbox.functions.proj_lineq(A=None, pinvA=None, **kwargs)

Bases: pyunlocbox.functions.proj

Projection on the plane satisfying the linear equality Az = y (eval, prox).

This function is the indicator function \(i_S(z)\) of the set \(S = \left\{z \in \mathbb{R}^N \mid Az = y \right\}\) that is zero if \(z\) is in the set and infinite otherwise.

The proximal operator is \(\operatorname{arg\,min}_z \| z - x \|_2 \text{ s.t. } Az = y\).

See generic attributes descriptions of the pyunlocbox.functions.proj base class. Note that the constructor takes keyword-only parameters.

See also

proj_b2

quadratic case

Notes

• A parameter pinvA, the pseudo-inverse of A, must be provided if the parameter A is provided as an operator/callable (not a matrix).

• The evaluation of this function is zero.

Examples

>>> from pyunlocbox import functions
>>> import numpy as np
>>> x = np.array([0, 0])
>>> A = np.array([[1, 1]])
>>> pinvA = np.linalg.pinv(A)
>>> y = np.array([1])
>>> f = functions.proj_lineq(A=A, pinvA=pinvA, y=y)
>>> sol = f.prox(x, 0)
>>> sol
array([0.5, 0.5])
>>> np.abs(A.dot(sol) - y) < 1e-15
array([ True])

class pyunlocbox.functions.structured_sparsity(lambda_=1, groups=[[]], weights=[0], **kwargs)

Bases: pyunlocbox.functions.func

Structured sparsity (eval, prox).

The structured sparsity term that is defined in the work of Jenatton et al. 2011 Proximal methods for hierarchical sparse coding.

\[\Omega(x) = \lambda \cdot \sum_{g \in G} w_g \cdot \|x_g\|_2\]

See generic attributes descriptions of the pyunlocbox.functions.func base class.

Parameters

lambda_float, optional

The scaling factor of the function that corresponds to \(\lambda\). Must be a non-negative number.

groups: list of lists of integers

Each element encodes the indices of the vector belonging to a single group. Corresponds to \(G\).

weightsarray_like

Weight associated to each group. Corresponds to \(w_g\). Must have the same length as \(G\).

Examples

>>> from pyunlocbox import functions
>>> groups = [[0, 1], [3, 2, 4]]
>>> weights = [2, 1]
>>> f = functions.structured_sparsity(10, groups, weights)
>>> x = [2, 2.5, -0.5, 0.3, 0.01]
>>> f.eval(x)
69.86305169905782
>>> f.prox(x, 0.1)
array([0.7506099 , 0.93826238, 0.        , 0.        , 0.        ])

Solvers

The pyunlocbox.solvers module implements a solving function (which will minimize your objective function) as well as common solvers.

Solving

Call solve() to solve your convex optimization problem using your instantiated solver and functions objects.

Interface

The solver base class defines a common interface to all solvers:

solver.pre(functions, x0)	Solver-specific pre-processing.
solver.algo(objective, niter)	Call the solver iterative algorithm and the provided acceleration scheme.
solver.post()	Solver-specific post-processing.

Solvers

Then, derived classes implement various common solvers.

gradient_descent(**kwargs)	Gradient descent algorithm.
forward_backward([accel])	Forward-backward proximal splitting (FISTA and ISTA) algorithm.
douglas_rachford([lambda_])	Douglas-Rachford proximal splitting algorithm.
generalized_forward_backward([lambda_])	Generalized forward-backward proximal splitting algorithm.

Primal-dual solvers (based on primal_dual)

mlfbf([L, Lt, d0])	Monotone+Lipschitz forward-backward-forward primal-dual algorithm.
projection_based([lambda_])	Projection-based primal-dual algorithm.

pyunlocbox.solvers.solve(functions, x0, solver=None, atol=None, dtol=None, rtol=0.001, xtol=None, maxit=200, verbosity='LOW', inplace=False)

Solve an optimization problem whose objective function is the sum of some convex functions.

This function minimizes the objective function \(f(x) = \sum\limits_{k=0}^{k=K} f_k(x)\), i.e. solves \(\operatorname{arg\,min}\limits_x f(x)\) for \(x \in \mathbb{R}^{n \times N}\) where \(n\) is the dimensionality of the data and \(N\) the number of independent problems. It returns a dictionary with the found solution and some informations about the algorithm execution.

Parameters

functionslist of objects

A list of convex functions to minimize. These are objects who must implement the pyunlocbox.functions.func.eval() method. The pyunlocbox.functions.func.grad() and / or pyunlocbox.functions.func.prox() methods are required by some solvers. Note also that some solvers can only handle two convex functions while others may handle more. Please refer to the documentation of the considered solver.

x0array_like

Starting point of the algorithm, \(x_0 \in \mathbb{R}^{n \times N}\).

solversolver class instance, optional

The solver algorithm. It is an object who must inherit from pyunlocbox.solvers.solver and implement the _pre(), _algo() and _post() methods. If no solver object are provided, a standard one will be chosen given the number of convex function objects and their implemented methods.

atolfloat, optional

The absolute tolerance stopping criterion. The algorithm stops when \(f(x^t) < atol\) where \(f(x^t)\) is the objective function at iteration \(t\). Default is None.

dtolfloat, optional

Stop when the objective function is stable enough, i.e. when \(\left|f(x^t) - f(x^{t-1})\right| < dtol\). Default is None.

rtolfloat, optional

The relative tolerance stopping criterion. The algorithm stops when \(\left|\frac{ f(x^t) - f(x^{t-1}) }{ f(x^t) }\right| < rtol\). Default is \(10^{-3}\).

xtolfloat, optional

Stop when the variable is stable enough, i.e. when \(\frac{\|x^t - x^{t-1}\|_2}{\sqrt{n N}} < xtol\). Note that additional memory will be used to store \(x^{t-1}\). Default is None.

maxitint, optional

The maximum number of iterations. Default is 200.

verbosity{‘NONE’, ‘LOW’, ‘HIGH’, ‘ALL’}, optional

The log level : 'NONE' for no log, 'LOW' for resume at convergence, 'HIGH' for info at all solving steps, 'ALL' for all possible outputs, including at each steps of the proximal operators computation. Default is 'LOW'.

inplacebool, optional

If True and x0 is a numpy array, then x0 will be modified in place during execution to save memory. It will then contain the solution. Be careful to pass data of the type (int, float32, float64) you want your computations to use.

Returns

solndarray

The problem solution.

solverstr

The used solver.

crit{‘ATOL’, ‘DTOL’, ‘RTOL’, ‘XTOL’, ‘MAXIT’}

The used stopping criterion. See above for definitions.

niterint

The number of iterations.

timefloat

The execution time in seconds.

objectivendarray

The successive evaluations of the objective function at each iteration.

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers

Define a problem:

>>> y = [4, 5, 6, 7]
>>> f = functions.norm_l2(y=y)

Solve it:

>>> x0 = np.zeros(len(y))
>>> ret = solvers.solve([f], x0, atol=1e-2, verbosity='ALL')
INFO: Dummy objective function added.
INFO: Selected solver: forward_backward
INFO: Forward-backward method
    dummy evaluation: 0.000000e+00
    norm_l2 evaluation: 1.260000e+02
Iteration 1 of forward_backward:
    dummy evaluation: 0.000000e+00
    norm_l2 evaluation: 1.400000e+01
    objective = 1.40e+01
Iteration 2 of forward_backward:
    dummy evaluation: 0.000000e+00
    norm_l2 evaluation: 2.963739e-01
    objective = 2.96e-01
Iteration 3 of forward_backward:
    dummy evaluation: 0.000000e+00
    norm_l2 evaluation: 7.902529e-02
    objective = 7.90e-02
Iteration 4 of forward_backward:
    dummy evaluation: 0.000000e+00
    norm_l2 evaluation: 5.752265e-02
    objective = 5.75e-02
Iteration 5 of forward_backward:
    dummy evaluation: 0.000000e+00
    norm_l2 evaluation: 5.142032e-03
    objective = 5.14e-03
Solution found after 5 iterations:
    objective function f(sol) = 5.142032e-03
    stopping criterion: ATOL

Verify the stopping criterion (should be smaller than atol=1e-2):

>>> np.linalg.norm(ret['sol'] - y)**2  
0.00514203...

Show the solution (should be close to y w.r.t. the L2-norm measure):

>>> ret['sol']
array([4.02555301, 5.03194126, 6.03832952, 7.04471777])

Show the used solver:

>>> ret['solver']
'forward_backward'

Show some information about the convergence:

>>> ret['crit']
'ATOL'
>>> ret['niter']
5
>>> ret['time']  
0.0012578964233398438
>>> ret['objective']  
[[126.0, 0], [13.99999999..., 0], [0.29637392..., 0], [0.07902528..., 0],
[0.05752265..., 0], [0.00514203..., 0]]

class pyunlocbox.solvers.solver(step=1.0, accel=None)

Bases: object

Defines the solver object interface.

This class defines the interface of a solver object intended to be passed to the pyunlocbox.solvers.solve() solving function. It is intended to be a base class for standard solvers which will implement the required methods. It can also be instantiated by user code and dynamically modified for rapid testing. This class also defines the generic attributes of all solver objects.

Parameters

stepfloat

The gradient-descent step-size. This parameter is bounded by 0 and \(\frac{2}{\beta}\) where \(\beta\) is the Lipschitz constant of the gradient of the smooth function (or a sum of smooth functions). Default is 1.

accelpyunlocbox.acceleration.accel

User-defined object used to adaptively change the current step size and solution while the algorithm is running. Default is a dummy object that returns unchanged values.

pre(functions, x0)

Solver-specific pre-processing. See parameters documentation in pyunlocbox.solvers.solve() documentation.

Notes

When preprocessing the functions, the solver should split them into two lists: * self.smooth_funs, for functions involved in gradient steps. * self.non_smooth_funs, for functions involved proximal steps. This way, any method that takes in the solver as argument, such as the methods in pyunlocbox.acceleration.accel, can have some context as to how the solver is using the functions.

algo(objective, niter)

Call the solver iterative algorithm and the provided acceleration scheme. See parameters documentation in pyunlocbox.solvers.solve()

Notes

The method self.accel.update_sol() is called before self._algo() because the acceleration schemes usually involves some sort of averaging of previous solutions, which can add some unwanted artifacts on the output solution. With this ordering, we guarantee that the output of solver.algo is not corrupted by the acceleration scheme.

Similarly, the method self.accel.update_step() is called after self._algo() to allow the step update procedure to act directly on the solution output by the underlying algorithm, and not on the intermediate solution output by the acceleration scheme in self.accel.update_sol().

post()

Solver-specific post-processing. Mainly used to delete references added during initialization so that the garbage collector can free the memory. See parameters documentation in pyunlocbox.solvers.solve().

objective(x)

Return the objective function at x.

Necessitate solver._pre(…) to be run first.

class pyunlocbox.solvers.gradient_descent(**kwargs)

Bases: pyunlocbox.solvers.solver

Gradient descent algorithm.

This algorithm solves optimization problems composed of the sum of any number of smooth functions.

See generic attributes descriptions of the pyunlocbox.solvers.solver base class.

Notes

This algorithm requires each function implement the pyunlocbox.functions.func.grad() method.

See pyunlocbox.acceleration.regularized_nonlinear for a very efficient acceleration scheme for this method.

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers
>>> dim = 25
>>> np.random.seed(0)
>>> xstar = np.random.rand(dim)  # True solution
>>> x0 = np.random.rand(dim)
>>> x0 = xstar + 5*(x0 - xstar) / np.linalg.norm(x0 - xstar)
>>> A = np.random.rand(dim, dim)
>>> step = 1 / np.linalg.norm(np.dot(A.T, A))
>>> f = functions.norm_l2(lambda_=0.5, A=A, y=np.dot(A, xstar))
>>> fd = functions.dummy()
>>> solver = solvers.gradient_descent(step=step)
>>> params = {'rtol':0, 'maxit':14000, 'verbosity':'NONE'}
>>> ret = solvers.solve([f, fd], x0, solver, **params)
>>> pctdiff = 100 * np.sum((xstar - ret['sol'])**2) / np.sum(xstar**2)
>>> print('Difference: {0:.1f}%'.format(pctdiff))
Difference: 1.3%

class pyunlocbox.solvers.forward_backward(accel=<pyunlocbox.acceleration.fista object>, **kwargs)

Bases: pyunlocbox.solvers.solver

Forward-backward proximal splitting (FISTA and ISTA) algorithm.

This algorithm solves convex optimization problems composed of the sum of a smooth and a non-smooth function.

See generic attributes descriptions of the pyunlocbox.solvers.solver base class.

Parameters

accelpyunlocbox.acceleration.accel

Acceleration scheme to use. Default is pyunlocbox.acceleration.fista(), which corresponds to the ‘FISTA’ solver. Passing pyunlocbox.acceleration.dummy() instead results in the ISTA solver. Note that while FISTA is much more time-efficient, it is less memory-efficient.

Notes

This algorithm requires one function to implement the pyunlocbox.functions.func.prox() method and the other one to implement the pyunlocbox.functions.func.grad() method.

See [BT09a] for details about the algorithm.

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers
>>> y = [4, 5, 6, 7]
>>> x0 = np.zeros(len(y))
>>> f1 = functions.norm_l2(y=y)
>>> f2 = functions.dummy()
>>> solver = solvers.forward_backward(step=0.5)
>>> ret = solvers.solve([f1, f2], x0, solver, atol=1e-5)
Solution found after 15 iterations:
    objective function f(sol) = 4.957288e-07
    stopping criterion: ATOL
>>> ret['sol']
array([4.0002509 , 5.00031362, 6.00037635, 7.00043907])

class pyunlocbox.solvers.generalized_forward_backward(lambda_=1, *args, **kwargs)

Bases: pyunlocbox.solvers.solver

Generalized forward-backward proximal splitting algorithm.

This algorithm solves convex optimization problems composed of the sum of any number of non-smooth (or smooth) functions.

See generic attributes descriptions of the pyunlocbox.solvers.solver base class.

Parameters

lambda_float, optional

A relaxation parameter bounded by 0 and 1. Default is 1.

Notes

This algorithm requires each function to either implement the pyunlocbox.functions.func.prox() method or the pyunlocbox.functions.func.grad() method.

See [RFP13] for details about the algorithm.

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers
>>> y = [0.01, 0.2, 8, 0.3, 0 , 0.03, 7]
>>> x0 = np.zeros(len(y))
>>> f1 = functions.norm_l2(y=y)
>>> f2 = functions.norm_l1()
>>> solver = solvers.generalized_forward_backward(lambda_=1, step=0.5)
>>> ret = solvers.solve([f1, f2], x0, solver)
Solution found after 2 iterations:
    objective function f(sol) = 1.463100e+01
    stopping criterion: RTOL
>>> ret['sol']
array([0. , 0. , 7.5, 0. , 0. , 0. , 6.5])

class pyunlocbox.solvers.douglas_rachford(lambda_=1, *args, **kwargs)

Bases: pyunlocbox.solvers.solver

Douglas-Rachford proximal splitting algorithm.

This algorithm solves convex optimization problems composed of the sum of two non-smooth (or smooth) functions.

See generic attributes descriptions of the pyunlocbox.solvers.solver base class.

Parameters

lambda_float, optional

The update term weight. It should be between 0 and 1. Default is 1.

Notes

This algorithm requires the two functions to implement the pyunlocbox.functions.func.prox() method.

See [CP07] for details about the algorithm.

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers
>>> y = [4, 5, 6, 7]
>>> x0 = np.zeros(len(y))
>>> f1 = functions.norm_l2(y=y)
>>> f2 = functions.dummy()
>>> solver = solvers.douglas_rachford(lambda_=1, step=1)
>>> ret = solvers.solve([f1, f2], x0, solver, atol=1e-5)
Solution found after 8 iterations:
    objective function f(sol) = 2.927052e-06
    stopping criterion: ATOL
>>> ret['sol']
array([3.99939034, 4.99923792, 5.99908551, 6.99893309])

class pyunlocbox.solvers.primal_dual(L=None, Lt=None, d0=None, *args, **kwargs)

Bases: pyunlocbox.solvers.solver

Parent class of all primal-dual algorithms.

See generic attributes descriptions of the pyunlocbox.solvers.solver base class.

Parameters

Lfunction or ndarray, optional

The transformation L that maps from the primal variable space to the dual variable space. Default is the identity, \(L(x)=x\). If L is an ndarray, it will be converted to the operator form.

Ltfunction or ndarray, optional

The adjoint operator. If Lt is an ndarray, it will be converted to the operator form. If L is an ndarray, default is the transpose of L. If L is a function, default is L, \(Lt(x)=L(x)\).

d0: ndarray, optional

Initialization of the dual variable.

class pyunlocbox.solvers.mlfbf(L=None, Lt=None, d0=None, *args, **kwargs)

Bases: pyunlocbox.solvers.primal_dual

Monotone+Lipschitz forward-backward-forward primal-dual algorithm.

This algorithm solves convex optimization problems with objective of the form \(f(x) + g(Lx) + h(x)\), where \(f\) and \(g\) are proper, convex, lower-semicontinuous functions with easy-to-compute proximity operators, and \(h\) has Lipschitz-continuous gradient with constant \(\beta\).

See generic attributes descriptions of the pyunlocbox.solvers.primal_dual base class.

Notes

The order of the functions matters: set \(f\) first on the list, \(g\) second, and \(h\) third.

This algorithm requires the first two functions to implement the pyunlocbox.functions.func.prox() method, and the third function to implement the pyunlocbox.functions.func.grad() method.

The step-size should be in the interval \(\left] 0, \frac{1}{\beta + \|L\|_{2}}\right[\).

See [KP15], Algorithm 6, for details.

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers
>>> y = np.array([294, 390, 361])
>>> L = np.array([[5, 9, 3], [7, 8, 5], [4, 4, 9], [0, 1, 7]])
>>> x0 = np.zeros(len(y))
>>> f = functions.dummy()
>>> f._prox = lambda x, T: np.maximum(np.zeros(len(x)), x)
>>> g = functions.norm_l2(lambda_=0.5)
>>> h = functions.norm_l2(y=y, lambda_=0.5)
>>> max_step = 1/(1 + np.linalg.norm(L, 2))
>>> solver = solvers.mlfbf(L=L, step=max_step/2.)
>>> ret = solvers.solve([f, g, h], x0, solver, maxit=1000, rtol=0)
Solution found after 1000 iterations:
    objective function f(sol) = 1.839060e+05
    stopping criterion: MAXIT
>>> ret['sol']
array([1., 1., 1.])

class pyunlocbox.solvers.projection_based(lambda_=1.0, *args, **kwargs)

Bases: pyunlocbox.solvers.primal_dual

Projection-based primal-dual algorithm.

This algorithm solves convex optimization problems with objective of the form \(f(x) + g(Lx)\), where \(f\) and \(g\) are proper, convex, lower-semicontinuous functions with easy-to-compute proximity operators.

See generic attributes descriptions of the pyunlocbox.solvers.primal_dual base class.

Parameters

lambda_float, optional

The update term weight. It should be between 0 and 2. Default is 1.

Notes

The order of the functions matters: set \(f\) first on the list, and \(g\) second.

This algorithm requires the two functions to implement the pyunlocbox.functions.func.prox() method.

The step-size should be in the interval \(\left] 0, \infty \right[\).

See [KP15], Algorithm 7, for details.

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers
>>> y = np.array([294, 390, 361])
>>> L = np.array([[5, 9, 3], [7, 8, 5], [4, 4, 9], [0, 1, 7]])
>>> x0 = np.array([500, 1000, -400])
>>> f = functions.norm_l1(y=y)
>>> g = functions.norm_l1()
>>> solver = solvers.projection_based(L=L, step=1.)
>>> ret = solvers.solve([f, g], x0, solver, maxit=1000, rtol=None, xtol=.1)
Solution found after 996 iterations:
    objective function f(sol) = 1.045000e+03
    stopping criterion: XTOL
>>> ret['sol']
array([0, 0, 0])

Acceleration

The pyunlocbox.acceleration module implements acceleration schemes for use with the pyunlocbox.solvers. Pass a given acceleration object as an argument to your chosen solver during its initialization so that the solver can use it.

Interface

The accel base class defines a common interface to all acceleration schemes:

accel.pre(functions, x0)	Pre-processing specific to the acceleration scheme.
accel.update_step(solver, objective, niter)	Update the step size for the next iteration.
accel.update_sol(solver, objective, niter)	Update the solution point for the next iteration.
accel.post()	Post-processing specific to the acceleration scheme.

Acceleration schemes

Then, derived classes implement various common acceleration schemes.

dummy()	Dummy acceleration scheme which does nothing.
backtracking([eta])	Backtracking acceleration based on a local quadratic approximation of the smooth part of the objective.
fista(**kwargs)	FISTA acceleration for forward-backward solvers.
fista_backtracking([eta])	FISTA acceleration with backtracking for forward-backward solvers.
regularized_nonlinear([k, lambda_, ...])	Regularized nonlinear acceleration (RNA) for gradient descent solvers.

class pyunlocbox.acceleration.accel

Bases: object

Defines the acceleration scheme object interface.

This class defines the interface of an acceleration scheme object intended to be passed to a solver inheriting from pyunlocbox.solvers.solver. It is intended to be a base class for standard acceleration schemes which will implement the required methods. It can also be instantiated by user code and dynamically modified for rapid testing. This class also defines the generic attributes of all acceleration scheme objects.

pre(functions, x0)

Pre-processing specific to the acceleration scheme.

Gets called when pyunlocbox.solvers.solve() starts running.

update_step(solver, objective, niter)

Update the step size for the next iteration.

Parameters

solverpyunlocbox.solvers.solver

Solver on which to act.

objectivelist of floats

List of evaluations of the objective function since the beginning of the iterative process.

niterint

Current iteration number.

Returns

float

Updated step size.

update_sol(solver, objective, niter)

Update the solution point for the next iteration.

Parameters

solverpyunlocbox.solvers.solver

Solver on which to act.

objectivelist of floats

List of evaluations of the objective function since the beginning of the iterative process.

niterint

Current iteration number.

Returns

array_like

Updated solution point.

post()

Post-processing specific to the acceleration scheme.

Mainly used to delete references added during initialization so that the garbage collector can free the memory. Gets called when pyunlocbox.solvers.solve() finishes running.

class pyunlocbox.acceleration.dummy

Bases: pyunlocbox.acceleration.accel

Dummy acceleration scheme which does nothing.

Used by default in most of the solvers. It simply returns unaltered the step size and solution point it receives.

class pyunlocbox.acceleration.backtracking(eta=0.5, **kwargs)

Bases: pyunlocbox.acceleration.dummy

Backtracking acceleration based on a local quadratic approximation of the smooth part of the objective.

Parameters

etafloat

A number between 0 and 1 representing the ratio of the geometric sequence formed by successive step sizes. In other words, it establishes the relation step_new = eta * step_old. Default is 0.5.

Notes

This is the backtracking strategy used in the original FISTA paper, [BT09a].

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers, acceleration
>>> y = [4, 5, 6, 7]
>>> x0 = np.zeros(len(y))
>>> f1 = functions.norm_l1(y=y, lambda_=1.0)
>>> f2 = functions.norm_l2(y=y, lambda_=0.8)
>>> accel = acceleration.backtracking()
>>> solver = solvers.forward_backward(accel=accel, step=10)
>>> ret = solvers.solve([f1, f2], x0, solver, atol=1e-32, rtol=None)
... 
Solution found after ... iterations:
    objective function f(sol) = 0.000000e+00
    stopping criterion: ATOL
>>> ret['sol']
array([4., 5., 6., 7.])

class pyunlocbox.acceleration.fista(**kwargs)

Bases: pyunlocbox.acceleration.dummy

FISTA acceleration for forward-backward solvers.

Notes

This is the acceleration scheme proposed in the original FISTA paper, [BT09a].

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers, acceleration
>>> y = [4, 5, 6, 7]
>>> x0 = np.zeros(len(y))
>>> f1 = functions.norm_l2(y=y)
>>> f2 = functions.dummy()
>>> accel = acceleration.fista()
>>> solver = solvers.forward_backward(accel=accel, step=0.5)
>>> ret = solvers.solve([f1, f2], x0, solver, atol=1e-5)
Solution found after 15 iterations:
    objective function f(sol) = 4.957288e-07
    stopping criterion: ATOL
>>> ret['sol']
array([4.0002509 , 5.00031362, 6.00037635, 7.00043907])

class pyunlocbox.acceleration.regularized_nonlinear(k=10, lambda_=1e-06, adaptive=True, dolinesearch=True, forcedecrease=True, **kwargs)

Bases: pyunlocbox.acceleration.dummy

Regularized nonlinear acceleration (RNA) for gradient descent solvers.

Parameters

kint, optional

Number of points to keep in the buffer for computing the extrapolation. (Default is 10.)

lambda_float or list of floats, optional

Regularization parameter in the acceleration scheme. The user can pass a list of candidates, and the acceleration algorithm will pick the one that provides the best extrapolation. (Default is 1e-6.)

adaptiveboolean, optional

If adaptive = True and the user has not provided a list of regularization parameters, the acceleration algorithm will assemble a grid of possible regularization parameters based on the SVD of the Gram matrix of vectors of differences in the extrapolation buffer. If adaptive = False, the algorithm will simply try to use the value(s) given in lambda_. (Default is True.)

dolinesearchboolean, optional

If dolinesearch = True, the acceleration scheme will try to return a point in the line segment between the current extrapolation and the previous one that provides a decrease in the value of the objective function. If dolinesearch = False, the algorithm simply returns the current extrapolation. (Default is True.)

forcedecreaseboolean, optional

If forcedecrese = True and we obtain a bad extrapolation, the algorithm returns the unchanged solution produced by the solver. If forcedecrease = False, the algorithm returns the new extrapolation no matter what. (Default is True.)

Notes

This is the acceleration scheme proposed in [SdB16].

See also Damien Scieur’s repository for the Matlab version that inspired this implementation.

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers, acceleration
>>> dim = 25
>>> np.random.seed(0)
>>> xstar = np.random.rand(dim)  # True solution.
>>> x0 = np.random.rand(dim)
>>> x0 = xstar + 5*(x0 - xstar) / np.linalg.norm(x0 - xstar)
>>> A = np.random.rand(dim, dim)
>>> step = 1 / np.linalg.norm(np.dot(A.T, A))
>>> f = functions.norm_l2(lambda_=0.5, A=A, y=np.dot(A, xstar))
>>> fd = functions.dummy()
>>> accel = acceleration.regularized_nonlinear(k=5)
>>> solver = solvers.gradient_descent(step=step, accel=accel)
>>> params = {'rtol':0, 'maxit':200, 'verbosity':'NONE'}
>>> ret = solvers.solve([f, fd], x0, solver, **params)
>>> pctdiff = 100 * np.sum((xstar - ret['sol'])**2) / np.sum(xstar**2)
>>> print('Difference: {0:.1f}%'.format(pctdiff))
Difference: 1.3%

property lambda_

class pyunlocbox.acceleration.fista_backtracking(eta=0.5, **kwargs)

Bases: pyunlocbox.acceleration.backtracking, pyunlocbox.acceleration.fista

FISTA acceleration with backtracking for forward-backward solvers.

Notes

This is the acceleration scheme and backtracking strategy proposed in the original FISTA paper, [BT09a].

Examples

>>> import numpy as np
>>> from pyunlocbox import functions, solvers, acceleration
>>> y = [4, 5, 6, 7]
>>> x0 = np.zeros(len(y))
>>> f1 = functions.norm_l2(y=y)
>>> f2 = functions.dummy()
>>> accel = acceleration.fista_backtracking()
>>> solver = solvers.forward_backward(accel=accel, step=0.5)
>>> ret = solvers.solve([f1, f2], x0, solver, atol=1e-5)
Solution found after 15 iterations:
    objective function f(sol) = 4.957288e-07
    stopping criterion: ATOL
>>> ret['sol']
array([4.0002509 , 5.00031362, 6.00037635, 7.00043907])

Operators

The pyunlocbox.operators module implements the following operators:

grad(x[, dim])	Returns the gradient of an array.
div(*args, **kwargs)	Returns the divergence of an array.

pyunlocbox.operators.grad(x, dim=2, **kwargs)

Returns the gradient of an array.

Parameters

dimint

Dimension of the grad

wxint

wyint

wzint

wtint

Weights to apply on each axis

Returns

dx, dy, dz, dtndarrays

Gradients following each axes, only the necessary ones are returned

Examples

>>> import numpy as np
>>> from pyunlocbox import operators
>>> x = np.arange(16).reshape(4, 4)
>>> dx, dy = operators.grad(x)

pyunlocbox.operators.div(*args, **kwargs)

Returns the divergence of an array.

Parameters

dxarray_like

dyarray_like

dzarray_like

dtarray_like

Arrays to operate on

Returns

xarray_like

Divergence vector

Examples

>>> import numpy as np
>>> from pyunlocbox import operators
>>> x = np.arange(16).reshape(4, 4)
>>> dx, dy = operators.grad(x)
>>> divx = operators.div(dx, dy)

Changelog

All notable changes to this project will be documented in this file. The format is based on Keep a Changelog and this project adheres to Semantic Versioning.

Unreleased

• New function: proj_lineq.

• New function: proj_sdsp.

• New function: proj_positive.

• New function: structured_sparsity.

• Continuous integration with Python 3.6, 3.7, 3.8, 3.9. Dropped 2.7, 3.4, 3.5.

• Merged all the extra requirements in a single dev requirement.

0.5.2 (2017-12-15)

Mostly a maintenance release. Much cleaning happened and a conda package is now available in conda-forge. Moreover, the package can now be tried online thanks to binder.

0.5.1 (2017-07-04)

Development status updated from Alpha to Beta.

New features:

• Acceleration module, decoupling acceleration strategies from the solvers

  • Backtracking scheme

  • FISTA acceleration

  • FISTA with backtracking

  • Regularized non-linear acceleration (RNA)

• Solvers: gradient descent algorithm

Bug fix:

• Decrease dimensionality of variables in Douglas Rachford tutorial to reduce test time and timeout on Travis CI.

Infrastructure:

• Continuous integration: dropped 3.3 (matplotlib dropped it), added 3.6

• We don’t build PDF documentation anymore. Less burden, HTML can be downloaded from readthedocs.

0.4.0 (2016-08-01)

New feature:

• Monotone+Lipschitz forward-backward-forward primal-dual algorithm (MLFBF)

Bug fix:

• Plots generated when building documentation (not stored in the repository)

Infrastructure:

• Continuous integration: dropped 2.6 and 3.2, added 3.5

• Travis-ci: check style and build doc

• Removed tox config (too cumbersome to use on dev box)

• Monitor code coverage and report to coveralls.io

0.3.0 (2015-05-29)

New features:

• Generalized forward-backward splitting algorithm

• Projection-based primal-dual algorithm

• TV-norm function (eval, prox)

• Nuclear-norm function (eval, prox)

• L2-norm proximal operator supports non-tight frames

• Two new tutorials using the TV-norm with Forward-Backward and Douglas-Rachford for image reconstruction and denoising

• New stopping criterion XTOL allows to stop when the variable is stable

Bug fix:

• Much more memory efficient. Note that the array which contains the initial solution is now modified in place.

0.2.1 (2014-08-20)

Bug fix version. Still experimental.

Bug fixes:

• Avoid complex casting to real

• Do not stop iterating if the objective function stays at zero

0.2.0 (2014-08-04)

Second usable version, available on GitHub and released on PyPI. Still experimental.

New features:

• Douglas-Rachford splitting algorithm

• Projection on the L2-ball for tight and non tight frames

• Compressed sensing tutorial using L2-ball, L2-norm and Douglas-Rachford

• Automatic solver selection

Infrastructure:

• Unit tests for all functions and solvers

• Continuous integration testing on Python 2.6, 2.7, 3.2, 3.3 and 3.4

0.1.0 (2014-06-08)

First usable version, available on GitHub and released on PyPI. Still experimental.

Features:

• Forward-backward splitting algorithm

• L1-norm function (eval and prox)

• L2-norm function (eval, grad and prox)

• Least square problem tutorial using L2-norm and forward-backward

• Compressed sensing tutorial using L1-norm, L2-norm and forward-backward

Infrastructure:

• Sphinx generated documentation using Numpy style docstrings

• Documentation hosted on Read the Docs

• Code hosted on GitHub

• Package hosted on PyPI

• Code checked by flake8

• Docstring and tutorial examples checked by doctest (as a test suite)

• Unit tests for functions module (as a test suite)

• All test suites executed in Python 2.6, 2.7 and 3.2 virtualenvs by tox

• Distributed automatic testing on Travis CI continuous integration platform

Contributing

Contributions are welcome, and they are greatly appreciated! The development of this package takes place on GitHub. Issues, bugs, and feature requests should be reported there. Code and documentation can be improved by submitting a pull request. Please add documentation and tests for any new code.

The package can be set up (ideally in a fresh virtual environment) for local development with the following:

$ git clone https://github.com/epfl-lts2/pyunlocbox.git
$ pip install --upgrade --editable pyunlocbox[dev]

The dev “extras requirement” ensures that dependencies required for development (to run the test suite and build the documentation) are installed.

You can improve or add solvers, functions, and acceleration schemes in pyunlocbox/solvers.py, pyunlocbox/functions.py, and pyunlocbox/acceleration.py, along with their corresponding unit tests in pyunlocbox/tests/test_*.py (with reasonable coverage). If you have a nice example to demonstrate the use of the introduced functionality, please consider adding a tutorial in doc/tutorials or a short example in examples.

Update README.rst and CHANGELOG.rst if applicable.

After making any change, please check the style, run the tests, and build the documentation with the following (enforced by Travis CI):

$ make lint
$ make test
$ make doc

Check the generated coverage report at htmlcov/index.html to make sure the tests reasonably cover the changes you’ve introduced.

To iterate faster, you can partially run the test suite, at various degrees of granularity, as follows:

$ python -m unittest pyunlocbox.tests.test_functions
$ python -m unittest pyunlocbox.tests.test_functions.TestCase.test_norm_l1

Making a release

1. Update the version number and release date in setup.py, pyunlocbox/__init__.py and CHANGELOG.rst.

2. Create a git tag with git tag -a v0.5.0 -m "PyUNLocBox v0.5.0".

3. Push the tag to GitHub with git push github v0.5.0. The tag should now appear in the releases and tags tab.

4. Create a release on GitHub and select the created tag. A DOI should then be issued by Zenodo.

5. Go on Zenodo and fix the metadata if necessary.

6. Build the distribution with make dist and check that the dist/pyunlocbox-0.5.0.tar.gz source archive contains all required files. The binary wheel should be found as dist/pyunlocbox-0.5.0-py2.py3-none-any.whl.

7. Test the upload and installation process:

   $ twine upload --repository-url https://test.pypi.org/legacy/ dist/*
   $ pip install --index-url https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple pyunlocbox
   

   Log in as the LTS2 user.

8. Build and upload the distribution to the real PyPI with make release.

9. Update the conda feedstock (at least the version number and sha256 in recipe/meta.yaml) by sending a PR to conda-forge.

References

BT09a

Amir Beck and Marc Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2009.

BT09b

Amir Beck and Marc Teboulle. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems. Image Processing, IEEE Transactions on, 2009.

CR07

Emmanuel Candes and Justin Romberg. Sparsity and incoherence in compressive sampling. Inverse problems, 2007. arXiv:math/0611957.

CP07

Patrick L Combettes and Jean-Christophe Pesquet. A douglas–rachford splitting approach to nonsmooth convex variational signal recovery. Selected Topics in Signal Processing, IEEE Journal of, 2007.

CP11

Patrick L Combettes and Jean-Christophe Pesquet. Proximal splitting methods in signal processing. In Fixed-Point Algorithms for Inverse Problems in Science and Engineering. 2011. arXiv:0912.3522.

KP15

Nikos Komodakis and Jean-Christophe Pesquet. Playing with duality. IEEE Signal Processing Magazine, 2015. arXiv:1406.5429.

RFP13

Hugo Raguet, Jalal Fadili, and Gabriel Peyré. A generalized forward-backward splitting. SIAM Journal on Imaging Sciences, 2013. arXiv:1108.4404.

SdB16

Damien Scieur, Alexandre dAspremont, and Francis Bach. Regularized nonlinear acceleration. arXiv, 2016. arXiv:1606.04133.