Cut-pursuit with preconditioned forward-Douglas-Rachford for regularization of classical functionals by graph total variation
Switch branches/tags
Nothing to show
Clone or download
1a7r0ch3
Latest commit 44b96df May 4, 2018

README.md

Cut-Pursuit with Preconditioned Forward-Douglas–Rachford for Minimizing Graph Total Variation with Additional Nondifferentiable Terms

Routines in C/C++.
Parallel implementation with OpenMP API.
MEX API for interface with GNU Octave or Matlab.
Boost API for interface with Python.

General problem statement

This extension of the cut-pursuit algorithm minimizes functionals structured over a weighted graph G = (V, E, w)

    F: xf(x) + ∑vV gv(xv) + ∑(u,v) ∈ E w(u,v)xuxv║ ,

where x ∈ ℍV for some base vector space ℍ, f is differentiable, and for all vV, gv admits directional derivatives on every points of its domain and every directions.

The cut-pursuit algorithm seeks partitions V of the set of vertices V, constituting the constant connected components of the solution, by successively solving the corresponding problem, structured over the reduced graph G = (V, E), that is

  arg minξ ∈ ℍV  F(x) ,    such that ∀ UV, ∀ uU, xu = ξU ,

and then refining the partition.
A key requirement is thus the ability to solve the reduced problem, which often have the exact same structure as the original one, but with much less vertices |V| ≪ |V|. If the solution of the original problem has only few constant connected components in comparison to the number of vertices, the cut-pursuit strategy can speed-up minimization by several orders of magnitude.
The preconditioned forward-Douglas–Rachford splitting algorithm is often well suited for such minimization, when the functionals are convex. We provide also MEX API for using it directly to solve the original problem; just remove the prefix CP_ in the methods listed below.

For the nonconvex case where the norm on the differences in the graph total variation is replaced by a ℓ0 norm, see L0 cut-pursuit repository, by Loïc Landrieu.

Available Routines

We provide implementations for a wide range of applications in convex cases, often used in signal processing or machine learning. We specify here the routines available with the MEX API, they are not all available in python yet. C/C++ routines are available separately (see below, § Files and Folders ).

Quadratic functional with ℓ1-norm regularization

The base space is ℍ = ℝ, and the general form is

    F: x ↦ 1/2 ║yAx2 + ∑vV λv |xv| + ∑(u,v) ∈ E w(u,v) |xuxv| ,

where y ∈ ℝn, and A: ℝn → ℝV is a linear operator, and λ ∈ ℝV and w ∈ ℝE are regularization weights.
This combination of ℓ1 norm and total variation is often coined fused LASSO.
Implementation also supports an additional positivity constraint   ∑vV ι+(xv), where ι+ is the convex indicator of ℝ+ : x ↦ 0 if x ≥ 0, +∞ otherwise.

Currently, A must be provided as a matrix.
There are several particular cases:

  • n ≪ |V|, call on CP_PFDR_graph_quadratic_d1_l1_mex
  • n ≥ |V|, call on CP_PFDR_graph_quadratic_d1_l1_AtA_mex with precomputations by A*
  • A is diagonal or f is a square weighted ℓ2 norm, call on CP_PFDR_graph_l22_d1_l1_mex

Quadratic functional with box constraints

The base space is ℍ = ℝ, and the general form is

    F: x ↦ 1/2 ║yA x2 + ∑vV ι[m,M](xv) + ∑(u,v) ∈ E w(u,v) |xuxv| ,

where y ∈ ℝn, and A: ℝn → ℝV is a linear operator, w ∈ ℝE are regularization weights, and ι[m,M] is the convex indicator of [m,M] : x ↦ 0 if mxM, +∞ otherwise.

Currently, A must be provided as a matrix.
There are several particular cases:

  • n ≪ |V|, call on CP_PFDR_graph_quadratic_d1_bounds_mex
  • n ≥ |V|, call on CP_PFDR_graph_quadratic_d1_bounds_AtA_mex with precomputations by A*
  • A is diagonal or f is a square weighted ℓ2 norm, call on CP_PFDR_graph_l22_d1_bounds_mex

Separable loss with simplex constraints

The base space is ℍ = ℝK, where K is a set of labels, and the general form is

    F: xf(y, x) + ∑vV ιΔK(xv) + ∑(u,v) ∈ E w(u,v)kK |xu,kxv,k| ,

where y ∈ ℝKV, f is a loss functional (see below), w ∈ ℝE are regularization weights, and ιΔK is the convex indicator of the simplex ΔK = {x ∈ ℝK | ∑k xk = 1 and ∀ k, xk ≥ 0}: x ↦ 0 if x ∈ ΔK, +∞ otherwise.

The following loss functionals are available, all implemented in the routine CP_PFDR_graph_loss_d1_simplex_mex:

  • linear, f(y, x) = − ∑vVkK xv,k yv,k
  • quadratic, f(y, x) = ∑vVkK (xv,kyv,k)2
  • smoothed Kullback–Leibler, f(y, x) = ∑vV KL(α u + (1 − α) yv, α u + (1 − α) xv),
    where α ∈ ]0,1[, u ∈ ΔK is the uniform discrete distribution over K, and KL: (p, q) ↦ ∑kK pk log(pk/qk).

Documentation

See § Available Routines for the problems currently implemented.

Directory tree

.   
├── data/       some illustrative data  
├── include/    C/C++ headers, with some doc  
├── octave/     GNU Octave or Matlab code  
│   ├── doc/    some documentation  
│   └── mex/    MEX API  
├── python/     python code and API  
└── src/        C/C++ sources  

C/C++

The C/C++ routines are documented within the corresponding headers in include/.

GNU Octave or Matlab

The MEX interfaces are documented within dedicated .m files in mex/doc/.
Beware that currently, inputs are not checked, wrong input types and sizes lead to segmentation faults or aberrant results.
See mex/compile_mex.m for typical compilation commands under UNIX systems.

The script example_EEG_CP.m exemplifies the problem § Quadratic functional with ℓ1-norm regularization, on a task of brain source identification with electroencephalography, described in the references.

Python

Currently, only the routine for the problem described in § Quadratic functional with ℓ1-norm regularization is implemented.
Some documentation can be found in the header of CP_quadratic_l1_py.cpp.

Make sure that your version of boost is at least 1.63, and that it is compiled with python support.
For compilation with cmake, the provided CMakeLists.txt file assumes python 3; this can be easily modified within the file. Compile with

cd python  
cmake .  
make   

You can compile the library in a given environment by replacing the cmake command by the following:

cmake . -DPYTHON_LIBRARY=$ENV/lib/libpython3.6m.so -DPYTHON_INCLUDE_DIR=$ENV/include/python3.6m -DBOOST_INCLUDEDIR=$ENV/include

with $ENV the path to the target environment (for anaconda it can be found with locate anaconda). Adapt the version of python to the one used in your environment (for anaconda you can find it with locate anaconda3/lib/libpython).

The script example_EEG_CP.py exemplifies the problem § Quadratic functional with ℓ1-norm regularization, on a task of brain source identification with electroencephalography, described in the references. The script requires the libraries matplotlib (tested with 2.2) and PyQt4.

ground truth raw retrieved activity identified sources

Data courtesy of Ahmad Karfoul and Isabelle Merlet, LTSI, INSERM U1099.

References

[1] H. Raguet and L. Landrieu, Cut-pursuit Algorithm for Regularizing Nonsmooth Functionals with Graph Total Variation.

[2] H. Raguet, A Note on the Forward-Douglas-Rachford Splitting Algorithm for Monotone Inclusion and Convex Optimization.