solvers: default debug prints for nonlinear solvers

[fnrizzi]

after refactoring the nonlinear solvers, default debug prints need to be re-added. need to decided which quantities to print

[fnrizzi]

@eparish1 I think we should also involve @pjb236 and @vbrunini and @jtencer to see what they would like to have printed

[pjb236]

Currently, we print FOM residual L2 norm, relative FOM residual L2 norm, ROM residual L2 norm, relative ROM residual L2 norm, and Gauss-Newton update L2 norm. I think we should keep all of these.

It might be nice to have some default printouts, as well as optional arguments for the ODE and/or solver objects to specify more printouts as needed. Perhaps this could take the form of some kind of "diagnostic" object/struct in which the user can create custom printouts. This would be useful for debugging and/or analysis. I am happy to expand on this idea if it sounds useful.

[fnrizzi]

@pjb236 "some kind of "diagnostic" object/struct in" yes this sounds a good idea! One thing to keep in mind is that for these prints, there should be no notion of FOM/ROM. These are just nonlinear LS solvers. So we should find the proper names

[pjb236]

@fnrizzi agreed on the names. We might want to look at other implementations of Gauss-Newton to determine good nomenclature, at least for the default print-outs. The user could specify whatever names they desire for the other norms in the diagnostic object. Here's a rough sketch of what I'm thinking:

`

struct Diagnostics{

// Constructor can initialize scaling for relative residuals, among other things. Diagnostics()

// operator method could be overloaded, templated as needed... void operator()(size_t step, scalar_t time, state_t x, residual_t r, ...) { // state_t should be the ROM state... // residual_t should be the full residual // perhaps we want other scalar inputs for norms that the solver already computes at each step.

// Users can compute norms, print data here as desired. }

} `

[fnrizzi]

yes, but that is a diagnostic for the ROM. We need to avoid any rom-related stuff for prints within the solvers as if ROMs are non-existing. We need a diagnostic scheme for solvers that is general.

[fnrizzi]

if you need ROM-related diagnostics, we can create anotehr isssue for that that must be addressed within the ROM package.

[pjb236]

No, everything I added is for a general nonlinear least squares system, I just called it the ROM state to illustrate it for our key problem of interest.

[fnrizzi]

I don't understand where the diagnostic object is passed and when. Here we are discussing debug print statements that are internal to pressio, not callbacks to classes passed-in. But we can also support that !

[pjb236]

Objects like this could be passed into the solver constructor, the ODE object constructor, or the ROM object constructor. It might be best to have different diagnostic objects for each of these, to allow for more fine grained control over what is printed. Also, this would be consistent with the modularity of Pressio. We would just have to make sure each diagnostic object does not interfere with the others.

[eparish1]

I have worked with a few CFD codes that have had this type of capability, and I think the cleanest implementation for this is in the form of a "hook". The hook gets called at, e.g., each time step, each nonlinear iteration. We come up with a list of reasonable inputs to the hook, and then the user can do whatever they want with them by writing their own hook. For example, printing out state norms, post processing. I believe that this should be done independently of the print statements, however.

[eparish1]

I would propose printing the following:

1.) Residual norm 2.) Rel. Residual norm 3.) Gradient norm. 4.) Rel. Gradient Norm 5.) Norm of correction 6.) Line search parameters (e.g., step size, regularization constant)

[fnrizzi]

@eparish1

• the prints you propose are for the default ones? this issue is just for default prints.
• please make sure you use the pressio prints commands with colors: #126
• the default prints should be added such that each stage (corrector, updater, looper) does self-contained prints without knowing what the others are printing.

[eparish1]

@fnrizzi "the default prints should be added such that each stage (corrector, updater, looper) does self-contained prints without knowing what the others are printing."

Can you explain more in detail what you mean by this? Are you thinking that each part will have separate print statements? Or alternatively each part could have separate print methods that are then brought together at the looper level.

[fnrizzi]

I am thinking (but have to see if this makes sense based on the code) that the looper would have a print method that prints information that it owns, and calls the print on the parent which prints its information and itself calls print on its parent (if it has one), etc... effectively, the mixins do not make sense as standalone pieces, so this chain is justifiable in my view. In other words, I don't think the bottom piece of the chain should print everything. every piece should only print the info it owns

[eparish1]

With this format, I think it will be hard to organize the print statements. For example, I think gradient and relative gradient should be printed next to each other. Relative gradient would be owned by the looper, and gradient would be owned by I believe the corrector. So we can't do something as coarse as

class looper{ void print(){ print looper owned statements corrector.print() }

We would have to do something more targeted like

class looper{ void print(){ print rel_grad corrector.print_grad() (and so on....) }

Am I making any sense, or am I missing your idea?

[fnrizzi]

ok, i see what you are saying... you are right!

I think I am getting convinced that we should have a separate free function that accepts a nonlinear solver object and is in charge of the prints. And this function is called from within the looper class. Something like this that lives somewhere:

template <typename solver_t, typename step_t>
void printNonlinearLeastSquaresMetrics(const solver_t & solver, step_t step)
{
  // get what you need from solver and print 
}

which we can always overload arbitrarily. The other advantage is that whatever changes we do to the solver classes, the print function does not need to change. Which is then callled as:

// for example, for the solve_until_max_iters:
template<typename system_t, typename state_t>
void solve(const system_t & sys, state_t & state)
{
  iteration_t iStep = 0;
  while (++iStep <= iterative_base_t::maxIters_)
  {
    T::computeCorrection(sys, state);
    T::updateState(sys, state);
#ifdef PRESSIO_ENABLE_DEBUG_PRINT
    printNonlinearLeastSquaresMetrics(*this, iStep);
#endif
  }
}

what do you think?

[eparish1]

I agree with this. I put things in as a void member of the looper class:

  void printSolverStatus(int iStep){
#ifdef PRESSIO_ENABLE_DEBUG_PRINT
      const auto & correction = T::viewCorrection();
      const auto correctionNorm = ::pressio::ops::norm2(correction);
      const auto & g = T::getGradient();
      const auto gnorm = ::pressio::ops::norm2(g);
      const auto resNorm = T::residualNormCurrentCorrectionStep();
      if (iStep == 1){
        gnorm0_ = ::pressio::ops::norm2(g) ;
        resNorm0_ = T::residualNormCurrentCorrectionStep();
      }

    ::pressio::utils::io::print_stdout(fmt, std::scientific,
            " Nonlinear iteration no. =" , iStep,
            utils::io::reset(),
            " Residual l2 norm (abs) =", resNorm,
            " Residual l2 norm (rel) =", resNorm/resNorm0_,
            " Gradient l2 norm (abs) =", gnorm,
            " Gradient l2 norm (rel) =", gnorm/gnorm0_,
            " Correction l2 norm =", correctionNorm,
            "\n");
#endif
  }

I do think that this would be better if it was it was a non-member function, or a class, as that we could reuse things across different solvers. So I'm doing this now.

[eparish1]

@fnrizzi What are your thoughts on this solution:

template <typename sc_t>
class NonlinearLeastSquaresDefaultMetricsPrinter{
private:
  sc_t norm0_ = {};

public:
  template <typename solver_t, typename step_t>
  void givenGradientNormsPrintRest(const solver_t & solver, step_t iStep, const sc_t & absGNorm, const sc_t & relGNorm)
  {   
      const auto correctionNorm = solver.correctionNorm();
      const auto resNorm = solver.residualNormCurrentCorrectionStep();
      if (iStep == 1) resNorm0_ = resNorm;
      printImpl(iStep, correctionNorm, resNorm, resNorm/norm0_, absGNorm, relGNorm);
  }

  template <typename solver_t, typename step_t>
  void givenResidualNormsPrintRest(const solver_t & solver, step_t iStep, const sc_t & absResNorm, const sc_t & relResNorm)
  {   
      const auto correctionNorm = solver.correctionNorm();
      const auto gNorm = solver.gradientNorm();
      if (iStep == 1) norm0_ = gNorm;
      printImpl(iStep, correctionNorm, absResNorm, relResNorm, gNorm, gNorm/norm0_);
  }

private:
  template <typename step_t>
  void printImpl(step_t iStep, 
                 const sc_t & correctionNorm,
                 const sc_t & absResNorm, const sc_t & relResNorm, 
                 const sc_t & absGNorm,   const sc_t & relGNorm) const 
  {
      auto fmt = utils::io::cyan() + utils::io::bold();      
      ::pressio::utils::io::print_stdout(fmt, std::scientific,
            " Nonlinear iteration no. =" , iStep,
            utils::io::reset(),
            " Residual l2 norm (abs) =", absResNorm,
            " Residual l2 norm (rel) =", relResNorm,
            " Gradient l2 norm (abs) =", absGNorm,
            " Gradient l2 norm (rel) =", relGNorm,
            " Correction l2 norm =", correctionNorm,
            "\n");
  }

};

template<bool absolute, typename sc_t, typename T>
class SolveUntilGradientNormBelowTol
  : public T,
    public IterativeBase< SolveUntilGradientNormBelowTol<absolute, sc_t, T>, sc_t >
{ 
  using this_t = SolveUntilGradientNormBelowTol<absolute, sc_t, T>;
  using iterative_base_t = IterativeBase<this_t, sc_t>;
  using typename iterative_base_t::iteration_t;

  iteration_t iStep_ = {};
  NonlinearLeastSquaresDefaultMetricsPrinter<sc_t> solverStatusPrinter = {};

public:
  SolveUntilGradientNormBelowTol() = delete;

  template <typename ...Args>
  SolveUntilGradientNormBelowTol(Args &&... args)
    : T(std::forward<Args>(args)...){}

  template<typename system_t, typename state_t>
  void solve(const system_t & sys, state_t & state)
  {
    sc_t initialNorm = {};
    sc_t absoluteNorm = {};
    sc_t relativeNorm = {};
    iStep_ = 0;
    while (++iStep_ <= iterative_base_t::maxIters_)
    {
      T::computeCorrection(sys, state);
      T::updateState(sys, state);

      const auto & g = T::getGradient();
      absoluteNorm = ::pressio::ops::norm2(g);
      if (iStep_==1) initialNorm = absoluteNorm;
      relativeNorm = absoluteNorm/absoluteNorm;

  #ifdef PRESSIO_ENABLE_DEBUG_PRINT
      solverStatusPrinter.givenGradientNormsPrintRest(*this, iStep_, absoluteNorm, relativeNorm);
  #endif 

      if (absolute){
        if (absoluteNorm < iterative_base_t::tolerance_)
          break;
      }
      else{
        if (relativeNorm < iterative_base_t::tolerance_)
          break;
      }
    }
  }

private:
  iteration_t getNumIterationsExecutedImpl() const {
    return iStep_;
  }

};

[fnrizzi]

@eparish1 I like the proposed version, thanks! I made a few edits:

• the printer should only be in charge of printing, not computing things. So we need to make sure the methods needed in the printer class (e.g. gradientNorm) are present in the nonlinear solver class

• i added two methods to the printer class, so that depending on what the looper class already has, it calls the appropriate one

• one things to check carefully and decide: who is in charge of holding/computing certain methods, e.g. gradientNorm. E.g., I think the norms of the operators (gradient, and residual) have to be computed no matter what. We need to decide who does that. I believe it should be the operator class at the top of the mixins chain, which already computes the residualNorm. Maybe we could do the same for the gradientNorm.

• We also need to be careful because the looper class does not know what method is used to compute the correction. For example, it could be the lin solve of the hessian/gradient system but also the QR method.

[eparish1]

@fnrizzi @pjb236 @jtencer When working with Python, they have a "verbosity" option, which controls how much information is printed. I think this would be nice. Thoughts?

[fnrizzi]

@eparish1

• I think for the time being, let's fix the default prints to be what we said above.
• but we can add also two more things later on:
  • support for accepting a callback for the print
  • your verbosity option: how is it defined in python? and is it defined for all solvers in python?

[eparish1]

@fnrizzi Sounds good. In python, verbosity usually takes on between a 1-4 input argument, 1 being the lease verbosity 4 being the most. It is implemented for most but not all solvers. The meaning of 1-4 varies across solvers.

[fnrizzi]

@eparish1 let's keep this issue just for default prints, I moved the verbosity to a new one #138

[eparish1]

@fnrizzi Per your comment on who owns the norms.

Norm of residual needs to be computed no matter what, but I don't think this is true for the gradient.

I'd just propose this, for example:

template<
  typename T, typename state_t, typename lin_solver_t, ::pressio::Norm normType
  >
class HessianGradientCorrector : public T
{
  using sc_t = typename ::pressio::containers::details::traits<state_t>::scalar_t;

  state_t correction_ = {};
  lin_solver_t & solverObj_;
  sc_t residualNorm_ = {};
  sc_t gradientNorm_ = {};
  sc_t correctionNorm_ = {};

public:
  static constexpr auto normType_ = normType;

  HessianGradientCorrector() = delete;

  template <typename system_t, typename ...Args>
  HessianGradientCorrector(const system_t & system,
         const state_t & state,
         lin_solver_t & solverObj,
         Args && ... args)
    : T(system, state, std::forward<Args>(args)...),
      correction_(state), solverObj_(solverObj){}

public:
  template <typename system_t>
  void computeCorrection(const system_t & sys, state_t & state)
  {
    T::computeOperators(sys, state, normType, residualNorm_);

    auto & H = T::getHessian();
    auto & g = T::getGradient();
    solverObj_.solveAllowMatOverwrite(H, g, correction_);

    correctionNorm_ = pressio::ops::norm2(correction_);
    gradientNorm_ = pressio::ops::norm2(g);
  }

  const state_t & viewCorrection() const{ return correction_; }
  const sc_t correctionNormFromCurrentCorrectionStep() const{ return correctionNorm_; }
  const sc_t gradientNormCurrentCorrectionStep() const{ return gradientNorm_; }
  const sc_t residualNormCurrentCorrectionStep() const{ return residualNorm_; }

  template< typename system_t>
  void residualNorm(const system_t & system, const state_t & state, sc_t & result){
    T::residualNorm(system, state, normType, result);
  }

};

[eparish1]

I pushed this to the branch issues/121. I am not sure how you want the printer class to be included. I put it in pressio_solvers.hpp.

[fnrizzi]

@eparish1 I like how it is shaping up!

• I agree on residual norm always needed, but not gradient norm. basically what you have now prints the gradient norm that is zero for the case that is not applicable, right? I think it is fine. We could also add something to indicate if a corrector supports gradient or not that we can detect from the printer class. but no big deal, maybe it is too complex for now.

• it makes sense that the corrector class has method to get the normsFromCurrentCorrectionStep. Just one small thing, make sure the names are consistent: blablaNormFromCurrentCorrectionStep, I think some methods currently lack the from.

• forgot one thing: let's also enclose the printer data member solverStatusPrinter inside the preproc directives, so that it is not even constructed if the PRESSIO_ENABLE_DEBUG_PRINT is off

• what do you think of the overall design we have come up with? do you like it? disagree with some things? anything you see we might be missing?

I pushed this to the branch issues/121. I am not sure how you want the printer class to be included. I put it in pressio_solvers.hpp.

Yes, that looks fine. The headers should be included such that the order makes sense because those global files pressio_ode.h` take care of including in the right order so that all headers do not need to have any include statements inside

[fnrizzi]

@eparish1

I think the current code:

T::computeCorrection(sys, state);
while(something)
{
  // print
  // check convergence 
  // update
  T::computeCorrection(sys, state);
}

iss the same as:

while(something)
{
  T::computeCorrection(sys, state);
  // print
  // check convergence 
  // update
}

right? Am i missing something? if you agree, i will switch to the second one since it is "cleaner" (at least to me).

Also, is there a reason you were getting references for the correction and gradient and recomputing norms inside the looping class? I changed ttht to get norms direcrly since they are already computed at the corrector step

[eparish1]

@fnrizzi

Good point! Yes this is the same :)

Per your second point, yes also good point. I think I did this before we decided to add these methods to the corrector step, and didn't realize the double computation. Thanks!