MIGRAD can call G2 even if HasG2 is set to false

[AdrianDBS]

Check duplicate issues.

• Checked for duplicates

Description

Minuit2's gradient-based minimizers are occasionally calling the hessian diagonal method G2 of a cost function even if the HasG2 property is set to false.
In other words, HasG2may just be ignored, at least during minimizer seeding.

See #20665

Reproducer

Example:

1. Create a FCNBase subclass for the Rosenbrock problem:

function = Pow(1 - x, 2) + 100*Pow(-Pow(x, 2) + y, 2);

gradient[x] = 400*x*(Pow(x, 2) - y) + 2*x - 2;
gradient[y] = -200*Pow(x, 2) + 200*y;

hessian[x, x] = 1200*Pow(x, 2) - 400*y + 2;
hessian[x, y] = -400*x;
hessian[y, x] = -400*x;
hessian[y, y] = 200;

which

• returns true for HasGradient and HasHessian
• returns the above-mentioned values for the cost function value, gradient and hessian
• returns false for HasG2
• throws an exception when G2 is called (or passively monitors a call to it)

2. Minimize that problem for the initial values (x=0, y=1) using MnMigrad
3. See G2 getting called although it shouldn't (at least to my mind)

ROOT version

ROOT v6.36.06

Installation method

Build from source

Operating system

Windows

Additional context

No response

[ferdymercury]

Thanks a lot for the report. Do you have a working reproducer (small ROOT script that we can run to reproduce) ?

[guitargeek]

No need to provide a reproducer. I have created one myself and I see the problem. I'm on it.

// root repro.C
#include <Minuit2/FCNBase.h>
#include <Minuit2/MnMigrad.h>
#include <Minuit2/MnUserParameters.h>
#include <Minuit2/FunctionMinimum.h>

#include <vector>
#include <iostream>
#include <cmath>

class RosenbrockFCN : public ROOT::Minuit2::FCNBase {

   static int gG2CallCount;

public:
   static int StaticG2CallCount() { return gG2CallCount; }

   double Up() const override { return 1.0; }

   // Rosenbrock: (1-x)^2 + 100*(y-x^2)^2
   double operator()(std::vector<double> const &p) const override
   {
      const double x = p.at(0);
      const double y = p.at(1);
      const double a = (1.0 - x);
      const double b = (y - x * x);
      return a * a + 100.0 * b * b;
   }

   bool HasGradient() const override { return true; }
   bool HasHessian() const override { return true; }

   // we claim we do NOT provide G2
   bool HasG2() const override { return false; }

   std::vector<double> Gradient(std::vector<double> const &p) const override
   {
      const double x = p.at(0);
      const double y = p.at(1);

      // gradient[x] = 400*x*(x^2 - y) + 2*x - 2
      // gradient[y] = -200*x^2 + 200*y
      std::vector<double> g(2);
      g[0] = 400.0 * x * (x * x - y) + 2.0 * x - 2.0;
      g[1] = -200.0 * x * x + 200.0 * y;
      return g;
   }

   // Hessian in packed upper-triangular storage (column-wise) like Minuit2 covariances:
   // for n=2: [H00, H01, H11]
   //
   // hxx = 1200*x^2 - 400*y + 2
   // hxy = -400*x
   // hyy = 200
   std::vector<double> Hessian(std::vector<double> const &p) const override
   {
      const double x = p.at(0);
      const double y = p.at(1);

      const double hxx = 1200.0 * x * x - 400.0 * y + 2.0;
      const double hxy = -400.0 * x;
      const double hyy = 200.0;

      std::vector<double> h(3);
      h[0] = hxx; // (0,0)
      h[1] = hxy; // (0,1)
      h[2] = hyy; // (1,1)
      return h;
   }

   std::vector<double> G2(std::vector<double> const &p) const override
   {
      ++gG2CallCount;
      std::vector<double> out(p.size());

      auto hessian = Hessian(p);
      out[0] = hessian[0];
      out[1] = hessian[2];

      return out;
   }
};

int RosenbrockFCN::gG2CallCount = 0;

void repro()
{
   RosenbrockFCN fcn;

   ROOT::Minuit2::MnUserParameters upar;
   // initial values and step size
   upar.Add("x", 0.0, 0.1);
   upar.Add("y", 1.0, 0.1);

   std::cout << "RosenbrockFCN::StaticG2CallCount() = " << fcn.StaticG2CallCount() << "\n";

   ROOT::Minuit2::MnMigrad migrad(fcn, upar);

   ROOT::Minuit2::FunctionMinimum min = migrad();

   std::cout << "Minimum valid: " << min.IsValid() << "\n";
   std::cout << "fval: " << min.Fval() << "\n";
   std::cout << "x = " << min.UserState().Value("x") << "\n";
   std::cout << "y = " << min.UserState().Value("y") << "\n";

   std::cout << "RosenbrockFCN::StaticG2CallCount() = " << fcn.StaticG2CallCount() << "\n";
}

The output:

RosenbrockFCN::StaticG2CallCount() = 0
Minimum valid: 1
fval: 0
x = 1
y = 1
RosenbrockFCN::StaticG2CallCount() = 1

[AdrianDBS]

@ferdymercury and @guitargeek, thanks for the fast response and for creating the reproducer. I would have otherwise tried to create one myself today.

[AdrianDBS]

@guitargeek, I just saw that you found a bug in the course of this ticket (#20936). I was just about to create another bug report that, I assume, will be fixed by your changes :-) I will double-check, though, once the fix is merged in.

[AdrianDBS]

One minor related question:
I saw that in your reproducer, you return only the upper-triangular elements for the Hessian. In other examples (e.g. in the test that was added for #20936), all elements (incl. the duplicate lower-triangular ones) are returned. Does the code auto-detect and accept both "output formats"?

[guitargeek]

No, you are right it needs to return the full n-by-n Hessian. I also realized that later