optimization, SolverBFGS

[Juan Carlos Araujo Cabarcas]

Dear all,

I would like to use SolverBFGS from deal.II/optimization/ in my project, and I try to follow the documentation in
    https://www.dealii.org/current/doxygen/deal.II/classSolverBFGS.html

First, I would like to mention that the documentation would really improve by adding a minimal example on how to use SolverBFGS.

In my project, I cannot just define an independent function to be pass to SolverBFGS.

Instead, I define a class ShapePolar, which needs to be initialized with some mesh parameters par, as in

    ShapePolar polar (par);
   
then I have implemented the function

    double ShapePolar::fun (const Vector<double>, Vector<double>)

that computes a scalar value of an objective function and its gradient Vector, when we pass X.


Then, my attempt, after looking at the documentation, is to do something like this:

  Vector<double> X (N);
     
  SolverControl                residual_control (N, 1e-7);
  //SolverBFGS::AdditionalData   data (10, true);

  SolverBFGS<Vector<double> > opt (residual_control);   //, data);

  ShapePolar  polar (par);
  std::function<double(const Vector<double>, Vector<double>)>  fun = polar.fun;

  opt.solve ( fun, X );

After several failed attempts (syntax errors), I realize that I am missing something fundamental here and could use some guidance.

I would appreciate any hint on how to achieve using SolverBFGS for the setting described above.

Kind regards,

    Juan C. Araújo Cabarcas

[Bruno Turcksin]

Juan,

Take a look at https://stackoverflow.com/questions/17131768/how-to-directly-bind-a-member-function-to-an-stdfunction-in-visual-studio-11

I advise using a lambda like the second reply suggests

Best,

Bruno

[Juan Carlos Araujo Cabarcas]

Thanks Bruno for your prompt reply,

I got the important part of the code compiling, however since I am new to this lambda (functional) functions, I do not know how to pass the new definitions to SolverBFGS.

Looking at your suggestion I do the following:

class Shape_compute {   
    public:   
        Parameters par;
   
        Shape_compute (Parameters epar):par(epar) {
            Register([=]( Vector<double> ex, Vector<double> dx ){ return objective_fun (ex,dx); });
        }
   
    void Register(std::function<double(Vector<double>, Vector<double>)> Callback) {} //fun = objective_fun;
   
    double objective_fun ( Vector<double> ex, Vector<double> dx) {
       
        vector<double> x (par.design_N), grad (par.design_N);
        for (unsigned int i=0; i < x.size (); i++) {
                x [i] = ex(i);
            }
        Adaptive::DtN_Helmholtz<2> fem (par, x, true);
            fem.run ();
            double objective_function = fem.objective_function;
           
            for (unsigned int i=0; i < par.design_N; i++) {
                if (i < par.design_N) {
                    vector<double> aux = x; aux [i] += h;
                    Adaptive::DtN_Helmholtz<2> fem (par, aux, false);
                    fem.run ();
                    grad [i] = (fem.objective_function - objective_function )/h;
                }
                dx(i) = grad [i];
            }
            iteration++;
            return objective_function;
        }
};

Then, when calling the optimization routine:

   
      Vector<double> X (N);
     
      SolverControl                residual_control (N, 1e-7);
      //SolverBFGS::AdditionalData   data (10, true);
     
      SolverBFGS<Vector<double> > opt (residual_control);//, data);
     

      Shape_compute  polar(par);
      opt.solve ( polar.Callback, X );     // If this line is commented, it compiles fine!

I would appreciate guidance on how to achieve this,

Thanks in advance,

  Juan Carlos Araújo

[Bruno Turcksin]

Here is a simple example:
https://wandbox.org/permlink/fbcldo9PIiHl6tfI In your case, try
something like

auto my_lambda = [&polar](Vector<double> a, Vector<double>b) {return
polar.Callback(a,b);};
opt.solve ( my_lambda, X );

Bruno

Le lun. 28 oct. 2019 à 06:30, Juan Carlos Araujo Cabarcas
<ju4...@gmail.com> a écrit :

[Juan Carlos Araujo Cabarcas]

Hi, thanks! that worked by using polar.objective_fun(a,b), instead of Callback.

Now, it seems that SolverBFGS is not starting, and terminates in the first function call!

However, the optimization should not be done since |grad x| is not small.

This is what I do:

     Vector<double> X (N), Dx (N);
     
      SolverControl                residual_control (100, 1e-10, true, true);
      residual_control.enable_history_data();
     
      SolverBFGS<Vector<double> > opt (residual_control,
                                                              SolverBFGS<Vector<double> >::AdditionalData(50, true) );
     
      Shape_compute  polar(par);

     
      auto my_lambda = [&polar](Vector<double> a, Vector<double>b) {

          return polar.objective_fun(a,b);
      };      //polar.Callback (a,b);};

            opt.solve ( my_lambda, X );
     

      std::cout << "          converged in " << residual_control.last_step()
          << " iterations" << std::endl;
     
      polar.objective_fun(X,Dx);     
      printf("\nThe norm of Dx is %.5f", Dx.l1_norm() );

and I get

          converged in 0 iterations

with The norm of Dx is 0.85597

[Bruno Turcksin]

There is a small bug in the lambda:

Le lun. 28 oct. 2019 à 09:36, Juan Carlos Araujo Cabarcas
<ju4...@gmail.com> a écrit :

> auto my_lambda = [&polar](Vector<double> a, Vector<double>b) {
> return polar.objective_fun(a,b);
> }; //polar.Callback (a,b);};

You want the vector to be passed by reference not by value, otherwise
the result cannot be updated. So it should be:
auto my_lambda = [&polar](const Vector<double> &a, Vector<double> &b) {
return polar.objective_fun(a,b);
};

Bruno

[Juan Carlos Araujo Cabarcas]

Thanks for the help, now it is running nicely!