iPiano
An implementation of the iPiano algorithms for non-convex and non-smooth optimization.
 All Classes Functions Variables
ipiano.h
1 
41 #ifndef IPIANO_H
42 #define IPIANO_H
43 
44 #include <random>
45 #include <fstream>
46 #include <functional>
47 #include <Eigen/Dense>
48 #include <glog/logging.h>
49 
50 #include "nmipiano.h"
51 
58 class iPiano : public nmiPiano
59 {
60 public:
61 
63  struct Options : public nmiPiano::Options
64  {
66  float beta_0m1 = 0.5f;
68  float c_1 = 1e-8;
70  float c_2 = 1e-8;
72  int steps = 10000;
73  };
74 
78  {
80  float beta_n;
82  float delta_n;
84  float gamma_n;
85  };
86 
90  static void silent_callback(const Iteration &iteration);
91 
95  static void default_callback(const Iteration &iteration, int n = 10);
96 
100  static void brief_callback(const Iteration &iteration, int n = 10);
101 
106  static void file_callback(const Iteration &iteration, const std::string &file);
107 
112  static void plot_callback(const Iteration &iteration, const std::string &file);
113 
120  iPiano(const std::function<float(const Eigen::MatrixXf&)> f,
121  const std::function<void(const Eigen::MatrixXf&, Eigen::MatrixXf&)> df,
122  std::function<float(const Eigen::MatrixXf&)> g,
123  const std::function<void(const Eigen::MatrixXf&, Eigen::MatrixXf&, float alpha)> prox_g,
124  Options options,
125  const std::function<void(const Iteration&)> callback = [](const Iteration &iteration) -> void { /* Nothing ... */ });
126 
128  ~iPiano();
129 
134  void optimize(Eigen::MatrixXf &x_star, float &f_x_star);
135 
136 protected:
137 
141  void initialize(Iteration &iteration);
142 
144  std::function<void(const Iteration&)> callback_;
147 };
148 
150 // iPiano::silent_callback
152 
154 {
155 
156 }
157 
159 // iPiano::default_callback
161 
162 void iPiano::default_callback(const iPiano::Iteration &iteration, int n)
163 {
164  if (iteration.n%n == 0)
165  {
166  std::cout << "[" << iteration.n << "] "
167  << iteration.f_x_n + iteration.g_x_n
168  << " (Delta_n = " << iteration.Delta_n << "; L_n = " << iteration.L_n
169  << "; alpha_n = " << iteration.alpha_n << "; beta_n = " << iteration.beta_n << ")"
170  << std::endl;
171  }
172 }
173 
175 // iPiano::default_callback
177 
178 void iPiano::brief_callback(const iPiano::Iteration &iteration, int n)
179 {
180  if (iteration.n%n == 0)
181  {
182  std::cout << "[" << iteration.n << "] "
183  << iteration.f_x_n + iteration.g_x_n
184  << " (Delta_n = " << iteration.Delta_n << ")"
185  << std::endl;
186  }
187 }
188 
190 // iPiano::file_callback
192 
193 void iPiano::file_callback(const iPiano::Iteration &iteration, const std::string &file)
194 {
195  std::ofstream out;
196  out.open(file, std::ios_base::app);
197  out << "iteration = " << iteration.n << "; f_x_n = " << iteration.f_x_n
198  << "; g_x_n = " << iteration.g_x_n << "; alpha_n = " << iteration.alpha_n
199  << "; beta_n = " << iteration.beta_n << "; L_n = " << iteration.L_n
200  << "; |df_x_n| = " << std::sqrt(iteration.df_x_n.squaredNorm())
201  << "; Delta_n = " << iteration.Delta_n << std::endl;
202  out.close();
203 }
204 
206 // iPiano::plot_callback
208 
209 void iPiano::plot_callback(const iPiano::Iteration &iteration, const std::string &file)
210 {
211  std::ofstream out;
212  out.open(file, std::ios_base::app);
213  out << iteration.n << "," << iteration.f_x_n + iteration.g_x_n
214  << "," << iteration.f_x_n << "," << iteration.g_x_n
215  << "," << iteration.alpha_n << "," << iteration.beta_n
216  << "," << iteration.L_n << "," << iteration.Delta_n << std::endl;
217  out.close();
218 }
219 
221 // iPiano::iPiano
223 
224 iPiano::iPiano(const std::function<float(const Eigen::MatrixXf&)> f,
225  const std::function<void(const Eigen::MatrixXf&,Eigen::MatrixXf&)> df,
226  std::function<float(const Eigen::MatrixXf &)> g,
227  const std::function<void(const Eigen::MatrixXf&, Eigen::MatrixXf&, float alpha)> prox_g,
228  iPiano::Options options,
229  const std::function<void(const iPiano::Iteration&)> callback)
230  : nmiPiano(f, df, g, prox_g, options)
231 {
232  options_ = options;
233  callback_ = callback;
234 }
235 
237 // iPiano::~iPiano
239 
241 {
242  // ...
243 }
244 
246 // iPiano::optimize
248 
249 void iPiano::optimize(Eigen::MatrixXf &x_star, float &f_x_star)
250 {
251  iPiano::Iteration iteration;
252  initialize(iteration);
253 
254  // Used for backtracking; referring to the next iterate x_np1:
255  Eigen::MatrixXf x_np1;
256 
257  for (unsigned int t = 0; t <= options_.max_iter; t++)
258  {
259  callback_(iteration);
260 
261  // Backtrack for the Lipschitz constant L_n and the updates
262  // alpha_n and beta_n:
263  float L_nm1 = estimateL(iteration);
264 
265  // Fall back to last L_n in worst case.
266  if (std::isinf(L_nm1) || std::isnan(L_nm1) || options_.BOUND_L_N)
267  {
268  LOG_IF(INFO, !options_.BOUND_L_N) << "Could not get starting value for local Lipschitz constant, using L_nm1 (L_n = " << L_nm1 << ")";
269  L_nm1 = iteration.L_n;
270  }
271 
272  float delta_nm1 = iteration.delta_n;
273  bool condition = false;
274 
275  int l = 0;
276  do
277  {
278  iteration.L_n = std::pow(options_.eta, l)*L_nm1;
279 
280  LOG_IF(FATAL, std::isinf(iteration.L_n) || std::isnan(iteration.L_n))
281  << "Could not find the local Lipschitz constant - L_n = "
282  << iteration.L_n << std::endl;
283  LOG_IF(INFO, l > 0 && l%1000 == 0)
284  << "Having a hard time finding the local Lipschitz constant (L_n = " << iteration.L_n
285  << "; L_nm1 = " << L_nm1 << "; eta = " << options_.eta << "; l = " << l << ")" << std::endl;
286 
287  float b = (delta_nm1 + iteration.L_n/2)/(options_.c_2 + iteration.L_n/2);
288  float beta = (b - 1)/(b - 1.f/2.f);
289 // float beta = 0.999;
290 
291  for (iteration.beta_n = beta; iteration.beta_n >= 0;
292  iteration.beta_n -= beta/options_.steps)
293  {
294  bool take = false; // Whether to take current beta_n.
295 // float alpha = 2*(1 - iteration.beta_n)/(iteration.L_n/2.f + options_.c_2);
296  float alpha = 2*(1 - iteration.beta_n)/(iteration.L_n);
297 
298  LOG_IF(FATAL, alpha < options_.c_1)
299  << "Cannot choose alpha_n - it does not exist: c_1 = "
300  << options_.c_1 << "; alpha_n = " << alpha
301  << " (L_n = " << iteration.L_n << ")!";
302 
303  for (iteration.alpha_n = alpha; iteration.alpha_n > options_.c_1;
304  iteration.alpha_n -= (alpha - options_.c_1)/options_.steps)
305  {
306  // TODO save some computation here!
307  iteration.delta_n = 1/iteration.alpha_n - iteration.L_n/2
308  - iteration.beta_n/(2*iteration.alpha_n);
309  iteration.gamma_n = 1/iteration.alpha_n - iteration.L_n/2
310  - iteration.beta_n/iteration.alpha_n;
311 
312  if (iteration.delta_n < delta_nm1
313  && iteration.delta_n >= iteration.gamma_n
314  && iteration.gamma_n >= options_.c_2)
315  {
316  // Good parameters, so take these!
317  take = true;
318  break;
319  }
320  }
321 
322  if (take)
323  {
324  break;
325  }
326  }
327 
328  // We fixes L_n and alpha_n, so we can try an iteration to check the
329  // Lipschitz condition.
330  iterate(iteration, iteration.beta_n, x_np1);
331  condition = checkLipschitz(iteration, x_np1);
332 
333  l++;
334  }
335  while (!condition); // Now alpha_n and L_n are set correctly.
336 
337  iteration.x_nm1 = iteration.x_n;
338  iteration.x_n = x_np1;
339 
340  // For iteration 0, this is done in initialize():
341  iteration.f_x_n = f_(iteration.x_n);
342  iteration.g_x_n = g_(iteration.x_n);
343  df_(iteration.x_n, iteration.df_x_n);
344 
345  LOG_IF(FATAL, iteration.x_n.rows() != iteration.df_x_n.rows()
346  || iteration.x_n.cols() != iteration.df_x_n.cols()) << "Output dimensions of df invalid.";
347 
348  iteration.n++;
349 
350  // Termination criterion
351  if (iteration.Delta_n < options_.epsilon && options_.epsilon > 0)
352  {
353  break;
354  }
355  }
356 
357  x_star = iteration.x_n;
358  f_x_star = iteration.f_x_n;
359 }
360 
362 // iPiano::initialize
364 
366 {
367  iteration.n = 0;
368  iteration.x_n = options_.x_0;
369  iteration.x_nm1 = options_.x_0;
370 
371  iteration.g_x_n = g_(iteration.x_n);
372  iteration.f_x_n = f_(iteration.x_n);
373 
374  df_(iteration.x_n, iteration.df_x_n);
375 
376  LOG_IF(FATAL, iteration.x_n.rows() != iteration.df_x_n.rows()
377  || iteration.x_n.cols() != iteration.df_x_n.cols()) << "Output dimensions of df invalid.";
378 
379  iteration.alpha_n = 0.1f; // Required for estimateL!
380 
381  // Estimate L_n.
382  if (options_.L_0m1 > 0)
383  {
384  iteration.L_n = options_.L_0m1;
385  }
386 
387  // Estimate L and take max.
388  if (iteration.df_x_n.squaredNorm() > 0.001)
389  {
390  iteration.L_n = std::max(iteration.L_n, estimateL(iteration));
391  }
392 
393  // beta_0m1 is given by the options!
394  iteration.beta_n = options_.beta_0m1;
395 
396  // Initialize alpha_0m1 to suffice alpha_0m1 <= 2*(1 - beta_0m)/L_0m1
397  // and delta_0m1 >= gamma_0m1 > c_2 where delta_0m1 and gamma_0m1 depend upon alpha_0m1.
398 
399 // float alpha = 2*(1 - iteration.beta_n)/(iteration.L_n/2.f + options_.c_2);
400  float alpha = 2*(1 - iteration.beta_n)/(iteration.L_n);
401 
402  LOG_IF(FATAL, alpha < options_.c_1)
403  << "Cannot choose alpha_n - it does not exist: c_1 = "
404  << options_.c_1 << "; alpha_n = " << alpha
405  << " (L_n = " << iteration.L_n << ")!";
406 
407  for (iteration.alpha_n = alpha; iteration.alpha_n > options_.c_1; iteration.alpha_n -= (alpha - options_.c_1)/options_.steps)
408  {
409  // TODO save some computation here!
410  iteration.delta_n = 1/iteration.alpha_n - iteration.L_n/2
411  - iteration.beta_n/(2*iteration.alpha_n);
412  iteration.gamma_n = 1/iteration.alpha_n - iteration.L_n/2
413  - iteration.beta_n/iteration.alpha_n;
414 
415  if (iteration.delta_n >= iteration.gamma_n && iteration.gamma_n >= options_.c_2)
416  {
417  // Good parameters, so take these!
418  break;
419  }
420  }
421 }
422 
423 #endif /* IPIANO_H */
424 
float L_n
L_n of current iterate (within iterations, this is also used to estimate L_np1 via backtracking)...
Definition: nmipiano.h:98
void initialize(Iteration &iteration)
Initialize the iteration structure (i.e. set iteration 0):
Definition: ipiano.h:365
std::function< float(const Eigen::MatrixXf &)> g_
Convex, nonsmooth part of objective.
Definition: nmipiano.h:182
float f_x_n
Smooth objective function at current iterate.
Definition: nmipiano.h:92
void optimize(Eigen::MatrixXf &x_star, float &f_x_star)
Optimize the given objective using the iPiano algorithm.
Definition: ipiano.h:249
Implementation of nmiPiano (algorithm 4) as proposed in [1]: [1] P. Ochs, Y. Chen, T. Brox, T. Pock. iPiano: Inertial Proximal Algorithm for Nonconvex Optimization. SIAM J. Imaging Sciences, vol. 7, no. 2, 2014.
Definition: nmipiano.h:56
int n
Current iteration.
Definition: nmipiano.h:102
Options of algorithm.
Definition: ipiano.h:63
iPiano(const std::function< float(const Eigen::MatrixXf &)> f, const std::function< void(const Eigen::MatrixXf &, Eigen::MatrixXf &)> df, std::function< float(const Eigen::MatrixXf &)> g, const std::function< void(const Eigen::MatrixXf &, Eigen::MatrixXf &, float alpha)> prox_g, Options options, const std::function< void(const Iteration &)> callback=[](const Iteration &iteration) -> void{})
Constructor.
Definition: ipiano.h:224
Structure representing an iteration, passed as is to a callback function to be able to monitor proces...
Definition: ipiano.h:77
float epsilon
Termination criterion; stop if Delta_n smaller than epsilon.
Definition: nmipiano.h:76
int steps
Number of dicsrete steps for alpha_n and beta_n to try.
Definition: ipiano.h:72
Eigen::MatrixXf df_x_n
Derivative of smooth objective at current iterate.
Definition: nmipiano.h:94
float c_2
Fixed c_2.
Definition: ipiano.h:70
float eta
Fixed eta for backtracking the local lipschitz constant.
Definition: nmipiano.h:70
float gamma_n
gamma_n of current iterate (within iterations, this is also used to estimate gamma_np1).
Definition: ipiano.h:84
std::function< void(const Eigen::MatrixXf &, Eigen::MatrixXf &)> df_
Gradient of smooth part of the objective function (i.e. derivative).
Definition: nmipiano.h:180
~iPiano()
Destructor.
Definition: ipiano.h:240
void iterate(Iteration &iteration, float beta, Eigen::MatrixXf &x_np1)
Do an iteration, i.e. compute x_np1.
Definition: nmipiano.h:425
float beta_n
beta_n of current iterate (within iterations, this is also used to estimate beta_np1).
Definition: ipiano.h:80
Eigen::MatrixXf x_nm1
Last iterate.
Definition: nmipiano.h:86
static void brief_callback(const Iteration &iteration, int n=10)
Brief callback to use with iPiano.
Definition: ipiano.h:178
float alpha_n
alpha_n of current iterate (within iterations, this is also used to estimate alpha_np1).
Definition: nmipiano.h:100
bool checkLipschitz(const Iteration &iteration, const Eigen::MatrixXf &x_np1)
Check the lipschitz condition for backtracking.
Definition: nmipiano.h:443
float g_x_n
Definition: nmipiano.h:96
Eigen::MatrixXf x_0
Initial iterate.
Definition: nmipiano.h:64
static void silent_callback(const Iteration &iteration)
Silent callback to use with iPiano.
Definition: ipiano.h:153
unsigned int max_iter
Maximum number of iterations.
Definition: nmipiano.h:66
static void plot_callback(const Iteration &iteration, const std::string &file)
File callback for plotting, writes CSV output of format: iteration,f,g,alpha_n,beta_n,L_n,Delta_n.
Definition: ipiano.h:209
float c_1
Fixed c_1.
Definition: ipiano.h:68
static void file_callback(const Iteration &iteration, const std::string &file)
File callback, used to write all information to the specified file.
Definition: ipiano.h:193
std::function< float(const Eigen::MatrixXf &)> f_
The smooth part of the objective function.
Definition: nmipiano.h:178
Implementation of iPiano (algorithm 5) as proposed in [1]: [1] P. Ochs, Y. Chen, T. Brox, T. Pock. iPiano: Inertial Proximal Algorithm for Nonconvex Optimization. SIAM J. Imaging Sciences, vol. 7, no. 2, 2014.
Definition: ipiano.h:58
std::function< void(const Iteration &)> callback_
Callback, note: overwrites nmiPiano::callback_!
Definition: ipiano.h:144
float Delta_n
Difference between two iterates.
Definition: nmipiano.h:90
Options of algorithm.
Definition: nmipiano.h:61
Eigen::MatrixXf x_n
Current iterate.
Definition: nmipiano.h:84
float delta_n
delta_n of current iterate (within iterations, this is also used to estimate delta_np1).
Definition: ipiano.h:82
float L_0m1
Initialization of loca Lipschitz.
Definition: nmipiano.h:72
float beta_0m1
beta_0m1 to initialize alpha_0m1, delta_0m1 and gamma_0m1 according to Equations (21) and (22)...
Definition: ipiano.h:66
bool BOUND_L_N
Whether to bound estimated Lipschitz constant below by the given L_n.
Definition: nmipiano.h:74
Options options_
Options (note: overwrites nmiPiano::options_!).
Definition: ipiano.h:146
Structure representing an iteration, passed as is to a callback function to be able to monitor proces...
Definition: nmipiano.h:81
static void default_callback(const Iteration &iteration, int n=10)
Default callback to use with iPiano.
Definition: ipiano.h:162