IPM alg nlp errors

Author: cnpetra

src/Optimization/InnerProduct.cpp

@@ -77,7 +77,7 @@ InnerProduct::~InnerProduct()
   delete vec_n_;
 }
 
-bool InnerProduct::apply_M(const hiopVector& x, hiopVector& y)
+bool InnerProduct::apply_M(const hiopVector& x, hiopVector& y) const
 {
   if(nlp_->useWeightedInnerProd()) {
     return nlp_->eval_M(x, y);
@@ -88,7 +88,7 @@ bool InnerProduct::apply_M(const hiopVector& x, hiopVector& y)
 }
 
 // Computes ||x||_M
-double InnerProduct::norm_M(const hiopVector& x)
+double InnerProduct::norm_M(const hiopVector& x) const
 {
   if(nlp_->useWeightedInnerProd()) {
     nlp_->eval_M(x, *vec_n_);
@@ -98,19 +98,19 @@ double InnerProduct::norm_M(const hiopVector& x)
     return x.twonorm();
   }
 }
-  // Computes H primal norm
-  double InnerProduct::norm_H_primal(const hiopVector& x)
-  {
-    if(nlp_->useWeightedInnerProd()) {      
-      nlp_->eval_H(x, *vec_n_);
-      auto dp = x.dotProductWith(*vec_n_);
+// Computes H primal norm
+double InnerProduct::norm_H_primal(const hiopVector& x) const
+{
+  if(nlp_->useWeightedInnerProd()) {      
+    nlp_->eval_H(x, *vec_n_);
+    auto dp = x.dotProductWith(*vec_n_);
       return ::std::sqrt(dp);
-    } else {
-      return x.twonorm();
-    }
+  } else {
+    return x.twonorm();
   }
+}
 // Computes H dual norm
-double InnerProduct::norm_H_dual(const hiopVector& x)
+double InnerProduct::norm_H_dual(const hiopVector& x) const
 {
   if(nlp_->useWeightedInnerProd()) {      
     nlp_->eval_H_inv(x, *vec_n_);
@@ -121,7 +121,7 @@ double InnerProduct::norm_H_dual(const hiopVector& x)
   }
 }
 
-double InnerProduct::norm_stationarity(const hiopVector& x)
+double InnerProduct::norm_stationarity(const hiopVector& x) const
 {
   if(nlp_->useWeightedInnerProd()) {
     nlp_->eval_H_inv(x, *vec_n_);
@@ -133,9 +133,10 @@ double InnerProduct::norm_stationarity(const hiopVector& x)
 }
 
 // Compute norm one weighted by M, i.e., 1^T*M*|x|
-double InnerProduct::norm_M_one(const hiopVector&x)
+double InnerProduct::norm_M_one(const hiopVector&x) const
 {
   if(nlp_->useWeightedInnerProd()) {
+    //opt! pre-compute M*1
     vec_n_->copyFrom(x);
     vec_n_->component_abs();
     nlp_->eval_M(*vec_n_, *vec_n2_);
@@ -146,4 +147,31 @@ double InnerProduct::norm_M_one(const hiopVector&x)
   }
 }
 
+double InnerProduct::norm_complementarity(const hiopVector& x) const
+{
+  if(nlp_->useWeightedInnerProd()) {
+    //opt! pre-compute M*1
+    vec_n_->copyFrom(x);
+    vec_n_->component_abs();
+    nlp_->eval_M(*vec_n_, *vec_n2_);
+    vec_n_->setToConstant(1.);
+    return vec_n_->dotProductWith(*vec_n2_);
+  } else {
+    return x.infnorm();
+  }
+}
+
+// Computes the "volume" of the space, 1^T M*1 
+double InnerProduct::volume() const
+{
+  if(nlp_->useWeightedInnerProd()) {
+    //opt! pre-compute M*1
+    vec_n_->setToConstant(1.);
+    nlp_->eval_M(*vec_n_, *vec_n2_);
+    return vec_n2_->dotProductWith(*vec_n_);
+  } else {
+    return nlp_->n_complem();
+  }
+}
+
 } // end namespace

src/Optimization/InnerProduct.hpp

@@ -87,30 +87,36 @@ class InnerProduct
   virtual ~InnerProduct();
 
   // Compute y=M*x
-  bool apply_M(const hiopVector& x, hiopVector& y);
+  bool apply_M(const hiopVector& x, hiopVector& y) const;
 
   // Computes ||x||_M
-  double norm_M(const hiopVector& x);
+  double norm_M(const hiopVector& x) const;
 
   // Computes H primal norm
-  double norm_H_primal(const hiopVector& x);
+  double norm_H_primal(const hiopVector& x) const;
 
   // Computes H dual norm
-  double norm_H_dual(const hiopVector& x);
+  double norm_H_dual(const hiopVector& x) const;
 
   // Computes norm of stationarity residual, using inf-norm for Euclidean spaces, H-inverse norm for Hilbert spaces
-  double norm_stationarity(const hiopVector& x);
+  double norm_stationarity(const hiopVector& x) const;
 
-  // Compute norm one weighted by M, i.e., 1^T*M*|x|
-  double norm_M_one(const hiopVector&x);
+  // Computes norm of complementarity, using inf-norm for Euclidean spaces, 1M norm for Hilbert spaces
+  double norm_complementarity(const hiopVector& x) const;
+  
+  // Computes 1M norm, i.e., ||x|| =  1^T*M*|x|
+  double norm_M_one(const hiopVector&x) const;
+
+  // Computes the "volume" of the space, norm of 1 function
+  double volume() const;
 private:
   // Pointer to "client" NLP
   hiopNlpFormulation* nlp_;
 
   // Working vector in the size n of the variables, allocated only when for the weighted case
-  hiopVector* vec_n_;
+  mutable hiopVector* vec_n_;
   // Working vector in the size n of the variables, allocated only when for the weighted case
-  hiopVector* vec_n2_;
+  mutable hiopVector* vec_n2_;
 };
 
 } //end namespace

src/Optimization/hiopAlgFilterIPM.cpp

@@ -656,10 +656,15 @@ bool hiopAlgFilterIPMBase::evalNlpAndLogErrors(const hiopIterate& it,
   // the one norms
   // double nrmDualBou=it.normOneOfBoundDuals();
   // double nrmDualEqu=it.normOneOfEqualityDuals();
-  double nrmDualBou, nrmDualEqu;
-  it.normOneOfDuals(nrmDualEqu, nrmDualBou);
+  double sc;
+  double sd;
+  double nrmDualBou;
+  double nrmDualEqu;
+  if(!it.compute_sc_sd(sc, sd, nrmDualEqu, nrmDualBou)) {
+    return false;
+  }
 
-  nlp->log->printf(hovScalars, "nrmOneDualEqu %g   nrmOneDualBo %g\n", nrmDualEqu, nrmDualBou);
+  nlp->log->printf(hovWarning, "nrmOneDualEqu %g   nrmOneDualBo %g\n", nrmDualEqu, nrmDualBou);
   if(nrmDualBou > 1e+10) {
     nlp->log->printf(hovWarning,
                      "Unusually large bound dual variables (norm1=%g) occured, "
@@ -684,8 +689,11 @@ bool hiopAlgFilterIPMBase::evalNlpAndLogErrors(const hiopIterate& it,
   }
 
   // scaling factors
-  double sd = fmax(p_smax, (nrmDualBou + nrmDualEqu) / (n + m)) / p_smax;
-  double sc = n == 0 ? 0 : fmax(p_smax, nrmDualBou / n) / p_smax;
+  //double sd = fmax(p_smax, (nrmDualBou + nrmDualEqu) / (n + m)) / p_smax;
+  //double sc = n == 0 ? 0 : fmax(p_smax, nrmDualBou / n) / p_smax;
+  
+  sd = fmax(p_smax, sd) / p_smax;
+  sc = n == 0 ? 0 : fmax(p_smax, sc) / p_smax;
 
   sd = fmin(sd, 1e+8);
   sc = fmin(sc, 1e+8);
@@ -696,7 +704,7 @@ bool hiopAlgFilterIPMBase::evalNlpAndLogErrors(const hiopIterate& it,
   // finally, the scaled nlp error
   nlpoverall = fmax(nlpoptim / sd, fmax(cons_violation, nlpcomplem / sc));
 
-  nlp->log->printf(hovScalars,
+  nlp->log->printf(hovWarning,
                    "nlpoverall %g  nloptim %g  sd %g  nlpfeas %g  nlpcomplem %g  sc %g cons_violation %g\n",
                    nlpoverall,
                    nlpoptim,
@@ -731,13 +739,15 @@ bool hiopAlgFilterIPMBase::evalNlpAndLogErrors2(const hiopIterate& it,
   nlp->runStats.tmSolverInternal.start();
 
   size_type n = nlp->n_complem(), m = nlp->m();
-  // the one norms
-  // double nrmDualBou=it.normOneOfBoundDuals();
-  // double nrmDualEqu=it.normOneOfEqualityDuals();
-  double nrmDualBou, nrmDualEqu;
-  it.normOneOfDuals(nrmDualEqu, nrmDualBou);
+  double sc;
+  double sd;
+  double nrmDualBou;
+  double nrmDualEqu;
+  if(!it.compute_sc_sd(sc, sd, nrmDualEqu, nrmDualBou)) {
+    return false;
+  }
 
-  nlp->log->printf(hovScalars, "nrmOneDualEqu %g   nrmOneDualBo %g\n", nrmDualEqu, nrmDualBou);
+  nlp->log->printf(hovWarning, "nrmOneDualEqu %g   nrmOneDualBo %g\n", nrmDualEqu, nrmDualBou);
   if(nrmDualBou > 1e+10) {
     nlp->log->printf(hovWarning,
                      "Unusually large bound dual variables (norm1=%g) occured, "
@@ -762,18 +772,21 @@ bool hiopAlgFilterIPMBase::evalNlpAndLogErrors2(const hiopIterate& it,
   }
 
   // scaling factors
-  //sd = max { p_smax, ||zl||_H + ||zu||_H + (||vl||_1 + ||vu||_1)/m } /p_smax
-  //sc = max { p_smax, ||zl||_H + ||zu||_H } /p_smax
+  //sd = max { p_smax, ||zl||_M + ||zu||_M + (||vl||_1 + ||vu||_1)/m } /p_smax
+  //sc = max { p_smax, ||zl||_M + ||zu||_M } /p_smax
   //nlpoptim = ||gradf + Jc'*yc + Jd'*Jd - M*zl - M*zu||_Hinv
   //           ||yd+vl-vu||_inf
   //nlpfeas = ||crhs- c||_inf
   //          ||drs- d||_inf
-  //          ||x-sxl-xl||_H
-  //          ||x-sxu+xu||_H
+  //          ||x-sxl-xl||_inf
+  //          ||x-sxu+xu||_inf
   //          ||d-sdl-dl||
   //          
-  double sd = fmax(p_smax, (nrmDualBou + nrmDualEqu) / (n + m)) / p_smax;
-  double sc = n == 0 ? 0 : fmax(p_smax, nrmDualBou / n) / p_smax;
+  //double sd = fmax(p_smax, (nrmDualBou + nrmDualEqu) / (n + m)) / p_smax;
+  //double sc = n == 0 ? 0 : fmax(p_smax, nrmDualBou / n) / p_smax;
+
+  sd = fmax(p_smax, sd) / p_smax;
+  sc = n == 0 ? 0 : fmax(p_smax, sc) / p_smax;
 
   sd = fmin(sd, 1e+8);
   sc = fmin(sc, 1e+8);
@@ -784,7 +797,7 @@ bool hiopAlgFilterIPMBase::evalNlpAndLogErrors2(const hiopIterate& it,
   // finally, the scaled nlp error
   nlpoverall = fmax(nlpoptim / sd, fmax(cons_violation, nlpcomplem / sc));
 
-  nlp->log->printf(hovScalars,
+  nlp->log->printf(hovWarning,
                    "nlpoverall %g  nloptim %g  sd %g  nlpfeas %g  nlpcomplem %g  sc %g cons_violation %g\n",
                    nlpoverall,
                    nlpoptim,
@@ -1159,7 +1172,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run()
   solver_status_ = NlpSolve_Pending;
 
   while(true) {
-    bret = evalNlpAndLogErrors(*it_curr,
+    bret = evalNlpAndLogErrors2(*it_curr,
                                *resid,
                                _mu,
                                _err_nlp_optim,
@@ -1252,7 +1265,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run()
 
       //! should perform only a partial update since NLP didn't change
       resid->update(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d, *logbar);
-      bret = evalNlpAndLogErrors(*it_curr,
+      bret = evalNlpAndLogErrors2(*it_curr,
                                  *resid,
                                  _mu,
                                  _err_nlp_optim,

src/Optimization/hiopIterate.cpp

@@ -236,26 +236,55 @@ double hiopIterate::normOneOfEqualityDuals() const
   return nrm1;
 }
 
-void hiopIterate::normOneOfDuals(double& nrm1Eq, double& nrm1Bnd) const
+// void hiopIterate::normOneOfDuals(double& nrm1Eq, double& nrm1Bnd) const
+// {
+// #ifdef HIOP_DEEPCHECKS
+//   assert(zl->matchesPattern(nlp->get_ixl()));
+//   assert(zu->matchesPattern(nlp->get_ixu()));
+//   assert(vl->matchesPattern(nlp->get_idl()));
+//   assert(vu->matchesPattern(nlp->get_idu()));
+// #endif
+//   // work locally with all the vectors. This will result in only one MPI_Allreduce call
+//   nrm1Bnd = zl->onenorm_local() + zu->onenorm_local();
+// #ifdef HIOP_USE_MPI
+//   double nrm1_global;
+//   int ierr = MPI_Allreduce(&nrm1Bnd, &nrm1_global, 1, MPI_DOUBLE, MPI_SUM, nlp->get_comm());
+//   assert(MPI_SUCCESS == ierr);
+//   nrm1Bnd = nrm1_global;
+// #endif
+//   nrm1Bnd += vl->onenorm_local() + vu->onenorm_local();
+//   nrm1Eq = yc->onenorm_local() + yd->onenorm_local();
+
+//   nlp->log->printf(hovWarning,
+//               "sc=--- sd=--- volume_n=--- nrm zl=%12.5e zu=%12.5e vl=%12.5e vu=%12.5e  yc=%12.5e yd=%12.5e\n",
+//                    zl->onenorm(), zu->onenorm(), vl->onenorm(), vu->onenorm(), yc->onenorm(), yd->onenorm());
+
+// }
+
+/// @brief Computes scaling factors sc and sd; also returns the "norms" of duals
+bool hiopIterate::compute_sc_sd(double& sc, double& sd, double& nrmDualsEq, double& nrmDualsBou) const
 {
-#ifdef HIOP_DEEPCHECKS
-  assert(zl->matchesPattern(nlp->get_ixl()));
-  assert(zu->matchesPattern(nlp->get_ixu()));
-  assert(vl->matchesPattern(nlp->get_idl()));
-  assert(vu->matchesPattern(nlp->get_idu()));
-#endif
-  // work locally with all the vectors. This will result in only one MPI_Allreduce call
-  nrm1Bnd = zl->onenorm_local() + zu->onenorm_local();
-#ifdef HIOP_USE_MPI
-  double nrm1_global;
-  int ierr = MPI_Allreduce(&nrm1Bnd, &nrm1_global, 1, MPI_DOUBLE, MPI_SUM, nlp->get_comm());
-  assert(MPI_SUCCESS == ierr);
-  nrm1Bnd = nrm1_global;
-#endif
-  nrm1Bnd += vl->onenorm_local() + vu->onenorm_local();
-  nrm1Eq = yc->onenorm_local() + yd->onenorm_local();
+  const auto nrmzl = nlp->inner_prod()->norm_M_one(*zl);
+  const auto nrmzu = nlp->inner_prod()->norm_M_one(*zu);
+  const auto volume_n = nlp->inner_prod()->volume();
+  const auto nrmvl = vl->onenorm_local();
+  const auto nrmvu = vu->onenorm_local();
+  nrmDualsBou = nrmzl + nrmzu + nrmvl + nrmvu;
+  sc = nrmDualsBou / volume_n;
+
+  const auto volume_m = nlp->m();
+  const auto nrmyc = yc->onenorm_local();
+  const auto nrmyd = yd->onenorm_local();
+  nrmDualsEq = nrmyc+nrmyd;
+  sd = (nrmDualsBou+nrmDualsEq)/(volume_n+volume_m);
+
+  nlp->log->printf(hovWarning,
+              "sc=%g sd=%g volume_n=%12.5e norms zl=%12.5e zu=%12.5e vl=%12.5e vu=%12.5e yc=%12.5e yd=%12.5e\n",
+                   sc, sd, volume_n, nrmzl, nrmzu, nrmvl, nrmvu, nrmyc, nrmyd);
+  return true;
 }
 
+  
 void hiopIterate::selectPattern()
 {
   sxl->selectPattern(nlp->get_ixl());

src/Optimization/hiopIterate.hpp

@@ -142,7 +142,10 @@ class hiopIterate
   virtual double normOneOfBoundDuals() const;
   virtual double normOneOfEqualityDuals() const;
   /* same as above but computed in one shot to save on communication and computation */
-  virtual void normOneOfDuals(double& nrm1Eq, double& nrm1Bnd) const;
+  //virtual void normOneOfDuals(double& nrm1Eq, double& nrm1Bnd) const;
+
+  /// @brief Computes scaling factors sc and sd; also returns the "norms" of duals
+  bool compute_sc_sd(double& sc, double& sd, double& nrmDualsEq, double& nrmDualsBou) const;
 
   /// @brief Entries corresponding to zeros in ix are set to zero
   virtual void selectPattern();

src/Optimization/hiopNlpFormulation.hpp

@@ -129,10 +129,11 @@ class hiopNlpFormulation
   /* Wrapper over user-defined inverse of H apply that also applies the NLP transformations. */
   bool eval_H_inv(const hiopVector& x, hiopVector& y); 
 
-  InnerProduct* inner_prod()
+  InnerProduct const* inner_prod() const
   {
     return inner_prod_;
   }
+
 protected:
   // calls specific hiopInterfaceXXX::eval_Jac_cons and deals with specializations of hiopMatrix arguments
   virtual bool eval_Jac_c_d_interface_impl(hiopVector& x, bool new_x, hiopMatrix& Jac_c, hiopMatrix& Jac_d) = 0;