iPiano
An implementation of the iPiano algorithms for nonconvex and nonsmooth optimization.

iPiano, proposed in [1], is an optimization algorithm combining forwardbackward splitting with an inertial force. This repository contains a C++ implementation of iPiano with applications to computer vision tasks. The implementation was submitted as part of a seminar paper [2] written at RWTH Aachen University and advised by Prof. Berkels.
[1] P. Ochs, Y. Chen, T. Brox, T. Pock. iPiano: Inertial Proximal Algorithm for Nonconvex Optimization SIAM Journal of Imaging Sciences, colume 7, number 2, 2014. [2] D. Stutz. Seminar paper "iPiano: Inertial Proximal Algorithm for NonConvex Optimization" https://github.com/davidstutz/seminaripiano
Similar to forwardbackward splitting, the algorithm tackles problems of the form
h(x) = f(x) + g(x)
where h and g are functions defined on \mathbb{R}^n
with different properties. In the most general setting, f is required to be smooth and g is required to be convex. The iterative algorithm is then described by the following update equation:
x^{(n + 1)} = \prox_{\alpha_n g}(x{(n)}  \nabla f(x{(n)}) + \beta_n (x^{(n)}  x^{(n  1)}))
where x^{(n)}
is the n
th iterate, \alpha_n
and \beta_n
are parameters, \prox_{\alpha_n g}
is the proximal mapping of g and \nabla f
is the gradient of f. Details can be found in [1] or the seminar paper corresponding to this implementation [2].
This implementation provides two variants of the algorithm: nmiPiano and iPiano. Ochs et al. [1] proved convergence of the latter one, while the former one is a simplified version. The algorithm can be applied to various tasks, as for example segmentation:
The project is based on CMake, Boost, Eigen, GLog as well as OpenCV (tested with OpenCV 2.x) and has been tested on Ubuntu 12.04 and Ubuntu 14.04:
sudo aptget install buildessential cmake libeigen3dev libboostalldev
For installation instructions of GLog and OpenCV, please see consult the respective web pages. The project is then compiled using:
git clone https://github.com/davidstutz/ipiano cd ipiano mkdir build cd build cmake .. make Scanning dependencies of target signal_denoising_cli [ 25%] Building CXX object signal_denoising_cli/CMakeFiles/signal_denoising_cli.dir/main.cpp.o Linking CXX executable signal_denoising_cli [ 25%] Built target signal_denoising_cli Scanning dependencies of target image_denoising_cli [ 50%] Building CXX object image_denoising_cli/CMakeFiles/image_denoising_cli.dir/main.cpp.o Linking CXX executable image_denoising_cli [ 50%] Built target image_denoising_cli Scanning dependencies of target phase_field_cli [ 75%] Building CXX object phase_field_cli/CMakeFiles/phase_field_cli.dir/main.cpp.o Linking CXX executable phase_field_cli [ 75%] Built target phase_field_cli Scanning dependencies of target phase_field_color_cli [100%] Building CXX object phase_field_color_cli/CMakeFiles/phase_field_color_cli.dir/main.cpp.o Linking CXX executable phase_field_color_cli [100%] Built target phase_field_color_cli
The modules in cmake/
may have to be adapted depending on the Eigen and GLog installations!
Usage can be illustrated using the example of onedimensional signal denoising as done in signal_denoising_cli
. The below code example shows the basic usage, also note the discussion below.
// We randomly sample a signal in the form of a Nx1 Eigen matrix: Eigen::MatrixXf signal; sampleSignal(signal); // perturbed_signal is the noisy signal we intend to denoise: Eigen::MatrixXf perturbed_signal = signal; perturbSignal(perturbed_signal); // Basic denoising functionals are provided in functionals.h // f_forentzianPairwise is a regularizer based on the lorentzian function; // i.e. it is differentiable but not convex. // g_absoluteUnary is an absolute data term which is convex but not differentiable. // The corresponding gradient and proximal mapping are implemented in functionals.h // and described in detail in [2]. std::function<float(const Eigen::MatrixXf&)> bound_f = std::bind(Functionals::Denoising::f_lorentzianPairwise, std::placeholders::_1, sigma, lambda); std::function<void(const Eigen::MatrixXf&, Eigen::MatrixXf&)> bound_df = std::bind(Functionals::Denoising::df_lorentzianPairwise, std::placeholders::_1, std::placeholders::_2, sigma, lambda); std::function<float(const Eigen::MatrixXf&)> bound_g = std::bind(Functionals::Denoising::g_absoluteUnary, std::placeholders::_1, perturbed_signal); std::function<void(const Eigen::MatrixXf&, Eigen::MatrixXf&, float)> bound_prox_g = std::bind(Functionals::Denoising::prox_g_absoluteUnary, std::placeholders::_1, perturbed_signal, std::placeholders::_2, std::placeholders::_3); // The initial iterate (i.e. starting point) will be a random signal. Eigen::MatrixXf x_0 = Eigen::MatrixXf::Zero(M, 1); perturbSignal(x_0); // nmiPiano provides the following options, see nmipiano.h or the discussion below. nmiPiano::Options nmi_options; nmi_options.x_0 = x_0; nmi_options.max_iter = 1000; nmi_options.L_0m1 = 100.f; nmi_options.beta = 0.5; nmi_options.eta = 1.05; nmi_options.epsilon = 1e8; // Both nmiPiano and iPiano provide callbacks to monitor progress. // The default_callback writes progress to std::cout, other callbacks to // write the progress to file are available in nmipiano.h and ipiano.h std::function<void(const nmiPiano::Iteration &iteration)> nmi_bound_callback = std::bind(nmiPiano::default_callback, std::placeholders::_1, 10); // For initialization, we provide f and g as well as their gradient/proximal mapping, // and the callback defined above. nmiPiano nmipiano(bound_f, bound_df, bound_g, bound_prox_g, nmi_options, nmi_bound_callback); // Optimization is done via .optimize expecting two arguments which will be the // final iterate as well as the corresponding function value (i.e. h = f + g) Eigen::MatrixXf nmi_x_star; float nmi_f_x_star; nmipiano.optimize(nmi_x_star, nmi_f_x_star);
The corresponding usage of iPiano can be found in signal_denoising_cli
. The functional to be optimized has to be provided as std::function
and is expected to have the following form:
static float f(const Eigen::MatrixXf &x); static void df_lorentzianPairwise(const Eigen::MatrixXf &x, Eigen::MatrixXf &df_x); static float g_absoluteUnary(const Eigen::MatrixXf &x); static void prox_g_absoluteUnary(const Eigen::MatrixXf &x, Eigen::MatrixXf &prox_f_x, float alpha);
Additional parameters are possible but have to be provided through std::bind
, as for example done with g_absoluteUnary
which in addition to the above parameters also expects the noisy signal in order to implement the absolute data term:
// Bind g_absoluteUnary such that the resulting std::function matches the above form! // Use std::placeholder::_1, std::placeholder::_2 etc. ... std::function<float(const Eigen::MatrixXf&)> bound_g = std::bind(Functionals::Denoising::g_absoluteUnary, std::placeholders::_1, perturbed_signal);
nmiPiano provides the following parameters (default values and documentation can also be found in nmipiano.h
):
x_0
: the initial iterate;max_iter
: maximum number of iterations;beta
: fixed \beta
, i.e. the parameter governing the momentum term/inertial force;eta
: parameter for backtracking to find the local Lipschitzconstant L_n in each iteration, see [1] or [2];L_0m1
: initial estimate of the local Lipschitzconstant; during initialization, the local Lipschitzconstant is estimated around x_0
and the maximum of the estimate and L_0m1
is taken;BOUND_L_N
: if true, the local Lipschitzconstant is always bounded below by L_0m1
;epsilon
: stopping criterion; if epsilon
is greater than zero, iterations stop if the squared norm of the difference of two consecutive iterates is smaller than epsilon
.Choosing these parameters needs some practice; reading [1] and/or [2] is highly recommended. In addition, the parameters strongly depend on the function to be optimized (e.g. if nmiPiano or iPiano do not converge for a given functional, try starting with a higher L_0m1
).
iPiano additionally provides the following parameters (also see ipiano.h
):
beta_0m1
: initial \beta
, i.e. the parameter governing the momentum term/inertial force  after the first iteration, \beta
is adapted automatically;c_1
: c_1
from [1] and [2], usually close to zero is fine, e.g. c_1 = 1e6
to c_1 = 1e12
;c_2
: same as c_1
;steps
: governs the resolution of finding appropriate \alpha_n
and \beta_n
in each iteration in order to guarantee convergence, see [2]; starting with high steps > 10000
is recommended  it can later be reduced depending on the functional.Independent of the function or the used variant, reading [1] and [2] is highly recommended!
Further functionals can be found in functionals.h
and another example can be found in image_denoising_cli
or phase_field_color_cli
.
This repository contains 4 examples for using nmiPiano and iPiano:
signal_denoising_cli
;image_denoising_cli
;phase_field_cli
;phase_field_color_cli
.The corresponding functions are detailed in [2]. Some examples are given below:
cd build cmake .. make # Opencv 4 windows containing the original, noisy and denoised signals! ./signal_denoising_cli/signal_denoising_cli [0] 79.8015 (Delta_n = 4.23812e38; L_n = 100; alpha_n = 0.01) [10] 50.5144 (Delta_n = 0.0272214; L_n = 653.026; alpha_n = 0.00153133) [20] 48.6092 (Delta_n = 0.0286435; L_n = 480.716; alpha_n = 0.00208023) [30] 46.641 (Delta_n = 0.0291762; L_n = 463.074; alpha_n = 0.00215948) [40] 44.6816 (Delta_n = 0.0290498; L_n = 461.657; alpha_n = 0.00216611) [50] 42.7609 (Delta_n = 0.0285756; L_n = 462.503; alpha_n = 0.00216215) # ... # Applies different functionals for denoising the image with added Gaussian noise: ./image_denoising_cli/image_denoising_cli ../3096.jpg [0] 16963.7 (Delta_n = 1.65681e37; L_n = 14.1274; alpha_n = 0.0707846) [10] 6109.42 (Delta_n = 1.39415; L_n = 209.762; alpha_n = 0.00476731) [20] 6004.02 (Delta_n = 0.078446; L_n = 180.711; alpha_n = 0.00553371) [30] 6002.3 (Delta_n = 0.0352952; L_n = 181.006; alpha_n = 0.00552467) [40] 6001.82 (Delta_n = 0.0147341; L_n = 181.044; alpha_n = 0.00552353) [50] 6001.69 (Delta_n = 0.0121218; L_n = 181.051; alpha_n = 0.00552331) # Applies a color phase field to segmentation (iteratively): ./phase_field_color_cli/phase_field_color_cli ../3096.jpg [0] 4.49466e+07 (Delta_n = 6.45593e38; L_n = 8.61429; alpha_n = 0.116086) [10] 4.48568e+07 (Delta_n = 34.6824; L_n = 8.81607; alpha_n = 0.113429) [20] 4.48357e+07 (Delta_n = 19.3647; L_n = 6.77839; alpha_n = 0.147528) [30] 4.48274e+07 (Delta_n = 13.7736; L_n = 7.1557; alpha_n = 0.139749) [40] 4.48222e+07 (Delta_n = 10.3676; L_n = 7.2935; alpha_n = 0.137108) [50] 4.48188e+07 (Delta_n = 9.12528; L_n = 6.54973; alpha_n = 0.152678) # ... [0] C_p = 0.43245,0.454655,0.534252; C_m = 0.44876,0.471163,0.554002 # ... [1] C_p = 0.342906,0.36354,0.445063; C_m = 0.488703,0.512195,0.593856 # ... [2] C_p = 0.303374,0.322877,0.399622; C_m = 0.479937,0.503389,0.586667 # ...
The provided example image is taken from the Berkeley Segmentation Dataset [3].
[3] P. Arbelaez, M. Maire, C. Fowlkes and J. Malik. Contour Detection and Hierarchical Image Segmentation. Transactions on Pattern Analysis and Machine Intelligence, volume 33, number 5, 2011.
Example segmentations are shown in the introduction.
Copyright (c) 2016, David Stutz All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.