Refinement on regularised system

Author: filikat

highs/ipm/hipo/factorhighs/Numeric.cpp

@@ -62,6 +62,11 @@ Int Numeric::solve(std::vector<double>& x) const {
   return kRetOk;
 }
 
-void Numeric::getReg(std::vector<double>& reg) { reg = std::move(total_reg_); }
+void Numeric::getReg(std::vector<double>& reg) {
+  // unpermute regularisation
+  permuteVector(total_reg_, S_->iperm());
+
+  reg = std::move(total_reg_);
+}
 
 }  // namespace hipo

highs/ipm/hipo/ipm/Iterate.cpp

@@ -627,52 +627,76 @@ void Iterate::finalResiduals(Info& info) const {
   }
 }
 
+void Iterate::setReg(LinearSolver& LS, OptionNla opt) {
+  // extract regularisation
+  LS.getReg(total_reg);
+
+  // easy access to primal/dual regularisation
+  if (opt == kOptionNlaNormEq) {
+    Rp = nullptr;
+    Rd = total_reg.data();
+  } else {
+    Rp = total_reg.data();
+    Rd = &total_reg[model->n()];
+  }
+}
+
 void Iterate::residuals6x6(const NewtonDir& d) {
   const std::vector<double>& dx = d.x;
   const std::vector<double>& dy = d.y;
   const std::vector<double>& dxl = d.xl;
   const std::vector<double>& dxu = d.xu;
   const std::vector<double>& dzl = d.zl;
   const std::vector<double>& dzu = d.zu;
+  const Int m = model->m();
+  const Int n = model->n();
+  assert(Rd);
 
   // res1,2,3,4,5,6 contain the rhs of the linear system
 
-  // ires1 = res1 - A * dx
+  // ires1 = res1 - A * dx - Rd * dy
   ires.r1 = res.r1;
   model->A().alphaProductPlusY(-1.0, dx, ires.r1);
+  for (Int i = 0; i < m; ++i) {
+    ires.r1[i] -= Rd[i] * dy[i];
+  }
 
   // ires2 = res2 - dx + dxl
-  for (Int i = 0; i < model->n(); ++i)
+  for (Int i = 0; i < n; ++i)
     if (model->hasLb(i))
       ires.r2[i] = res.r2[i] - dx[i] + dxl[i];
     else
       ires.r2[i] = 0.0;
 
   // ires3 = res3 - dx - dxu
-  for (Int i = 0; i < model->n(); ++i)
+  for (Int i = 0; i < n; ++i)
     if (model->hasUb(i))
       ires.r3[i] = res.r3[i] - dx[i] - dxu[i];
     else
       ires.r3[i] = 0.0;
 
-  // ires4 = res4 - A^T * dy - dzl + dzu
+  // ires4 = res4 - A^T * dy - dzl + dzu + Rp * dx
   ires.r4 = res.r4;
-  for (Int i = 0; i < model->n(); ++i) {
+  for (Int i = 0; i < n; ++i) {
     if (model->hasLb(i)) ires.r4[i] -= dzl[i];
     if (model->hasUb(i)) ires.r4[i] += dzu[i];
   }
   model->A().alphaProductPlusY(-1.0, dy, ires.r4, true);
+  for (Int i = 0; i < n; ++i) {
+    double reg_p = Rp ? Rp[i] : regul.primal;
+    ires.r4[i] += reg_p * dx[i];
+  }
 
   // ires5 = res5 - zl * dxl - xl * dzl
-  for (Int i = 0; i < model->n(); ++i) {
+  for (Int i = 0; i < n; ++i) {
     if (model->hasLb(i))
       ires.r5[i] = res.r5[i] - zl[i] * dxl[i] - xl[i] * dzl[i];
     else
       ires.r5[i] = 0.0;
   }
 
   // ires6 = res6 - zu * dxu - xu * dzu
-  for (Int i = 0; i < model->n(); ++i) {
+  for (Int i = 0; i < n; ++i) {
     if (model->hasUb(i))
       ires.r6[i] = res.r6[i] - zu[i] * dxu[i] - xu[i] * dzu[i];
     else

highs/ipm/hipo/ipm/Iterate.h

@@ -5,6 +5,7 @@
 
 #include "Info.h"
 #include "IpmData.h"
+#include "LinearSolver.h"
 #include "Model.h"
 #include "ipm/hipo/auxiliary/IntConfig.h"
 
@@ -55,6 +56,9 @@ struct Iterate {
 
   // regularisation values
   Regularisation& regul;
+  std::vector<double> total_reg;
+  double* Rp;
+  double* Rd;
 
   // residuals of linear system for iterative refinement
   Residuals ires;
@@ -191,6 +195,8 @@ struct Iterate {
   // Compute residual of 6x6 linear system for iterative refinement.
   // ===================================================================================
   void residuals6x6(const NewtonDir& d);
+
+  void setReg(LinearSolver& LS, OptionNla opt);
 };
 
 }  // namespace hipo

highs/ipm/hipo/ipm/LinearSolver.h

@@ -26,6 +26,7 @@ namespace hipo {
 // - flops: return number of flops needed for factorisation
 // - spops: return number of sparse ops needed for factorisation
 // - nz: return number of nonzeros in factorisation
+// - getReg: extract regularisation
 //
 // NB: forming the normal equations or augmented system is delegated to the
 // linear solver chosen, so that only the appropriate data (upper triangle,

highs/ipm/hipo/ipm/Refine.cpp

@@ -81,19 +81,26 @@ double Solver::computeOmega(const NewtonDir& delta) const {
 
   // tau_i =
   // tau_const * (inf_norm_rows(big 6x6 matrix)_i * inf_norm_delta + |rhs_i|)
+
   const double tau_const = 1000.0 * (5 * n_ + m_) * 1e-16;
   double omega_1{}, omega_2{};
 
+  assert(it_->Rd);
+
   // First block
   for (Int i = 0; i < m_; ++i) {
-    const double tau = tau_const * (model_.infNormRows(i) * inf_norm_delta +
-                                    std::abs(it_->res.r1[i]));
-    updateOmega(tau, omega_1, omega_2,
-                // numerator
-                std::abs(it_->ires.r1[i]),
-                // denominators
-                abs_prod_A[i], std::abs(it_->res.r1[i]),
-                model_.oneNormRows(i) * inf_norm_delta);
+    const double tau =
+        tau_const * (std::max(model_.infNormRows(i), std::abs(it_->Rd[i])) *
+                         inf_norm_delta +
+                     std::abs(it_->res.r1[i]));
+    updateOmega(
+        tau, omega_1, omega_2,
+        // numerator
+        std::abs(it_->ires.r1[i]),
+        // denominators
+        abs_prod_A[i] + std::abs(it_->Rd[i]) * std::abs(delta.y[i]),
+        std::abs(it_->res.r1[i]),
+        (model_.oneNormRows(i) + std::abs(it_->Rd[i])) * inf_norm_delta);
   }
 
   // Second block
@@ -128,20 +135,25 @@ double Solver::computeOmega(const NewtonDir& delta) const {
 
   // Fourth block
   for (Int i = 0; i < n_; ++i) {
+    double reg_p = it_->Rp ? it_->Rp[i] : regul_.primal;
+
     const double tau =
-        tau_const * (std::max(model_.infNormCols(i), 1.0) * inf_norm_delta +
-                     std::abs(it_->res.r4[i]));
+        tau_const *
+        (std::max(std::max(model_.infNormCols(i), 1.0), std::abs(reg_p)) *
+             inf_norm_delta +
+         std::abs(it_->res.r4[i]));
 
     double denom1 = abs_prod_At[i];
     if (model_.hasLb(i)) denom1 += std::abs(delta.zl[i]);
     if (model_.hasUb(i)) denom1 += std::abs(delta.zu[i]);
+    denom1 += std::abs(reg_p) * std::abs(delta.x[i]);
 
     updateOmega(tau, omega_1, omega_2,
                 // numerator
                 std::abs(it_->ires.r4[i]),
                 // denominators
                 denom1, std::abs(it_->res.r4[i]),
-                (model_.oneNormCols(i) + 2) * inf_norm_delta);
+                (model_.oneNormCols(i) + 2 + std::abs(reg_p)) * inf_norm_delta);
   }
 
   // Fifth block

highs/ipm/hipo/ipm/Solver.cpp

@@ -293,12 +293,14 @@ bool Solver::solve2x2(NewtonDir& delta, const Residuals& rhs) {
     std::vector<double> res8 = it_->residual8(rhs, res7);
 
     // factorise normal equations, if not yet done
-    if (!LS_->valid_)
+    if (!LS_->valid_) {
       if (Int status = LS_->factorNE(model_.A(), theta_inv)) {
         logH_.printe("Error while factorising normal equations\n");
         info_.status = (Status)status;
         return true;
       }
+      it_->setReg(*LS_, options_.nla);
+    }
 
     // solve with normal equations
     if (Int status = LS_->solveNE(res8, delta.y)) {
@@ -321,12 +323,14 @@ bool Solver::solve2x2(NewtonDir& delta, const Residuals& rhs) {
   // AUGMENTED SYSTEM
   else {
     // factorise augmented system, if not yet done
-    if (!LS_->valid_)
+    if (!LS_->valid_) {
       if (Int status = LS_->factorAS(model_.A(), theta_inv)) {
         logH_.printe("Error while factorising augmented system\n");
         info_.status = (Status)status;
         return true;
       }
+      it_->setReg(*LS_, options_.nla);
+    }
 
     // solve with augmented system
     if (Int status = LS_->solveAS(res7, rhs.r1, delta.x, delta.y)) {