Merge pull request #1669 from ERGO-Code/quass_steepest_edge

Quass: Steepest edge pricing, common reinvert

cmake/sources-python.cmake

@@ -369,6 +369,7 @@ set(highs_headers_python
     src/qpsolver/devexpricing.hpp
     src/qpsolver/eventhandler.hpp
     src/qpsolver/factor.hpp
+    src/qpsolver/feasibility_bounded.hpp
     src/qpsolver/feasibility_highs.hpp
     src/qpsolver/feasibility_quass.hpp
     src/qpsolver/feasibility.hpp

cmake/sources.cmake

@@ -372,6 +372,7 @@ set(highs_headers
     qpsolver/devexpricing.hpp
     qpsolver/eventhandler.hpp
     qpsolver/factor.hpp
+    qpsolver/feasibility_bounded.hpp
     qpsolver/feasibility_highs.hpp
     qpsolver/feasibility_quass.hpp
     qpsolver/feasibility.hpp

src/lp_data/Highs.cpp

@@ -3462,6 +3462,20 @@ HighsStatus Highs::callSolveQp() {
   settings.iterationlimit = options_.simplex_iteration_limit;
   settings.lambda_zero_threshold = options_.dual_feasibility_tolerance;
 
+  switch (options_.simplex_primal_edge_weight_strategy) {
+    case 0:
+      settings.pricing = PricingStrategy::DantzigWolfe;
+      break;
+    case 1:
+      settings.pricing = PricingStrategy::Devex;
+      break;
+    case 2:
+      settings.pricing = PricingStrategy::SteepestEdge;
+      break;
+    default:
+      settings.pricing = PricingStrategy::Devex;
+  }
+
   // print header for QP solver output
   highsLogUser(options_.log_options, HighsLogType::kInfo,
                "Iteration, Runtime, ObjVal, NullspaceDim\n");

src/qpsolver/a_asm.cpp

@@ -2,8 +2,7 @@
 #include "qpsolver/quass.hpp"
 
 QpAsmStatus solveqp_actual(Instance& instance, Settings& settings, QpHotstartInformation& startinfo, Statistics& stats, QpModelStatus& status, QpSolution& solution, HighsTimer& qp_timer) {
-  Runtime rt(instance);
-  rt.statistics = stats;
+  Runtime rt(instance, stats);
   rt.settings = settings;
   Quass quass(rt);
 

src/qpsolver/a_quass.cpp

@@ -2,6 +2,7 @@
 #include "qpsolver/a_asm.hpp"
 
 #include "qpsolver/feasibility_highs.hpp"
+#include "qpsolver/feasibility_bounded.hpp"
 
 
 
@@ -26,10 +27,21 @@ QpAsmStatus solveqp(Instance& instance, Settings& settings, Statistics& stats, Q
 
   // compute initial feasible point
   QpHotstartInformation startinfo(instance.num_var, instance.num_con);
-  computestartingpoint_highs(instance, settings, stats, modelstatus, startinfo, qp_timer);
-  if (modelstatus == QpModelStatus::INFEASIBLE) {
-    return QpAsmStatus::OK;
-  }
+  if (instance.num_con == 0 && instance.num_var <= 15000) {
+    computestartingpoint_bounded(instance, settings, stats, modelstatus, startinfo, qp_timer);
+    if (modelstatus == QpModelStatus::OPTIMAL) {
+      solution.primal = startinfo.primal;
+      return QpAsmStatus::OK;
+    }
+    if (modelstatus == QpModelStatus::UNBOUNDED) {
+      return QpAsmStatus::OK;
+    }
+  } else  {
+    computestartingpoint_highs(instance, settings, stats, modelstatus, startinfo, qp_timer);
+    if (modelstatus == QpModelStatus::INFEASIBLE) {
+      return QpAsmStatus::OK;
+    }
+  } 
 
   // solve
   QpAsmStatus status = solveqp_actual(instance, settings, startinfo, stats, modelstatus, solution, qp_timer);

src/qpsolver/basis.cpp

@@ -88,6 +88,7 @@ void Basis::rebuild() {
        i < activeconstraintidx.size() + nonactiveconstraintsidx.size(); i++) {
     constraintindexinbasisfactor[baseindex[i]] = i;
   }
+  reinversion_hint = false;
 }
 
 void Basis::report() {
@@ -166,7 +167,7 @@ void Basis::updatebasis(const Settings& settings, HighsInt newactivecon, HighsIn
 
   updatessinceinvert++;
   if (updatessinceinvert >= settings.reinvertfrequency || hint != 99999) {
-    rebuild();
+    reinversion_hint = true;
   }
   // since basis changed, buffered values are no longer valid
   buffered_p = -1;

src/qpsolver/basis.hpp

@@ -73,7 +73,6 @@ class Basis {
   std::vector<HighsInt> constraintindexinbasisfactor;
 
   void build();
-  void rebuild();
 
   // buffer to avoid recreating vectors
   Vector buffer_column_aq;
@@ -85,10 +84,17 @@ class Basis {
   HVector row_ep;
   HVector col_aq;
 
+  bool reinversion_hint = false;
  public:
+
+
   Basis(Runtime& rt, std::vector<HighsInt> active,
         std::vector<BasisStatus> atlower, std::vector<HighsInt> inactive);
 
+  bool getreinversionhint() { return reinversion_hint; }
+
+  void rebuild();
+
   HighsInt getnupdatessinceinvert() { return updatessinceinvert; }
 
   HighsInt getnumactive() const { return activeconstraintidx.size(); };

src/qpsolver/dantzigpricing.hpp

@@ -59,6 +59,10 @@ class DantzigPricing : public Pricing {
     return minidx;
   }
 
+  void recompute() {
+    // do nothing
+  }
+
   void update_weights(const Vector& aq, const Vector& ep, HighsInt p,
                       HighsInt q) {
     // does nothing

src/qpsolver/devexharrispricing.hpp

@@ -63,6 +63,10 @@ class DevexHarrisPricing : public Pricing {
     return minidx;
   }
 
+  void recompute() {
+    // do nothing
+  }
+
   void update_weights(const Vector& aq, const Vector& ep, HighsInt p,
                       HighsInt q) {
     HighsInt rowindex_p = basis.getindexinfactor()[p];

src/qpsolver/devexpricing.hpp

@@ -73,6 +73,10 @@ class DevexPricing : public Pricing {
     return minidx;
   }
 
+  void recompute() {
+    // do nothing
+  }
+
   void update_weights(const Vector& aq, const Vector& ep, HighsInt p,
                       HighsInt q) {
     HighsInt rowindex_p = basis.getindexinfactor()[p];

src/qpsolver/factor.hpp

@@ -26,6 +26,33 @@ class CholeskyFactor {
   bool has_negative_eigenvalue = false;
   std::vector<double> a;
 
+  void resize(HighsInt new_k_max) {
+    std::vector<double> L_old = L;
+    L.clear();
+    L.resize((new_k_max) * (new_k_max));
+    const HighsInt l_size = L.size();
+    // Driven by #958, changes made in following lines to avoid array
+    // bound error when new_k_max < current_k_max
+    HighsInt min_k_max = min(new_k_max, current_k_max);
+    for (HighsInt i = 0; i < min_k_max; i++) {
+      for (HighsInt j = 0; j < min_k_max; j++) {
+	assert(i * (new_k_max) + j < l_size);
+        L[i * (new_k_max) + j] = L_old[i * current_k_max + j];
+      }
+    }
+    current_k_max = new_k_max;
+  }
+
+ public:
+  CholeskyFactor(Runtime& rt, Basis& bas) : runtime(rt), basis(bas) {
+    uptodate = false;
+    current_k_max =
+        max(min((HighsInt)ceil(rt.instance.num_var / 16.0), (HighsInt)1000),
+            basis.getnuminactive());
+    L.resize(current_k_max * current_k_max);
+  }
+
+
   void recompute() {
     std::vector<std::vector<double>> orig;
     HighsInt dim_ns = basis.getinactive().size();
@@ -70,32 +97,6 @@ class CholeskyFactor {
     uptodate = true;
   }
 
-  void resize(HighsInt new_k_max) {
-    std::vector<double> L_old = L;
-    L.clear();
-    L.resize((new_k_max) * (new_k_max));
-    const HighsInt l_size = L.size();
-    // Driven by #958, changes made in following lines to avoid array
-    // bound error when new_k_max < current_k_max
-    HighsInt min_k_max = min(new_k_max, current_k_max);
-    for (HighsInt i = 0; i < min_k_max; i++) {
-      for (HighsInt j = 0; j < min_k_max; j++) {
-	assert(i * (new_k_max) + j < l_size);
-        L[i * (new_k_max) + j] = L_old[i * current_k_max + j];
-      }
-    }
-    current_k_max = new_k_max;
-  }
-
- public:
-  CholeskyFactor(Runtime& rt, Basis& bas) : runtime(rt), basis(bas) {
-    uptodate = false;
-    current_k_max =
-        max(min((HighsInt)ceil(rt.instance.num_var / 16.0), (HighsInt)1000),
-            basis.getnuminactive());
-    L.resize(current_k_max * current_k_max);
-  }
-
   QpSolverStatus expand(const Vector& yp, Vector& gyp, Vector& l, Vector& m) {
     if (!uptodate) {
       return QpSolverStatus::OK;
@@ -203,7 +204,7 @@ class CholeskyFactor {
   }
 
   void solve(Vector& rhs) {
-    if (!uptodate || (numberofreduces >= runtime.instance.num_con / 2 &&
+    if (!uptodate || (numberofreduces >= runtime.instance.num_var / 2 &&
                       !has_negative_eigenvalue)) {
       recompute();
     }

src/qpsolver/feasibility.hpp

@@ -3,6 +3,7 @@
 
 #include "runtime.hpp"
 #include "crashsolution.hpp"
+#include "feasibility_bounded.hpp"
 #include "feasibility_highs.hpp"
 #include "feasibility_quass.hpp"
 
@@ -16,6 +17,9 @@ void computestartingpoint(Runtime& runtime, QpHotstartInformation& result) {
         case Phase1Strategy::QUASS:
             computestartingpoint_quass(runtime, result);
             break;
+        case Phase1Strategy::BOUNDED:
+            computestartingpoint_bounded(runtime.instance, runtime.settings, runtime.statistics, runtime.status, result, qp_timer);
+            break;
     }
 }
 

src/qpsolver/feasibility_bounded.hpp

@@ -0,0 +1,100 @@
+#ifndef __SRC_LIB_FEASIBILITYBOUNDED_HPP__
+#define __SRC_LIB_FEASIBILITYBOUNDED_HPP__
+
+#include "Highs.h"
+#include "qpsolver/a_asm.hpp"
+#include "qpsolver/crashsolution.hpp"
+
+static void computestartingpoint_bounded(Instance& instance, Settings& settings, Statistics& stats, QpModelStatus& modelstatus, QpHotstartInformation& result, HighsTimer& timer) {
+  // compute initial feasible point for problems with bounds only (no general linear constraints)
+
+  // compute  Qx + c = 0 --> x = Q^-1c
+  std::vector<double> L;
+  L.resize(instance.num_var * instance.num_var);
+
+  // compute cholesky factorization of Q
+  for (size_t col = 0; col < instance.num_var; col++) {
+    for (size_t idx = instance.Q.mat.start[col]; idx < instance.Q.mat.start[col+1]; idx++) {
+      double sum = 0;
+      size_t row = instance.Q.mat.index[idx]; 
+      if (row == col) {
+        for (size_t k = 0; k < row; k++)
+          sum += L[k * instance.num_var + row] * L[k * instance.num_var + row];
+        L[row * instance.num_var + row] = sqrt(instance.Q.mat.value[idx] - sum);
+      } else {
+        for (size_t k = 0; k < row; k++)
+          sum += (L[k * instance.num_var + col] * L[k * instance.num_var + row]);
+        L[row * instance.num_var + col] =
+            (instance.Q.mat.value[idx] - sum) / L[row * instance.num_var + row];
+      }
+    }
+  }
+
+  // solve for c
+  Vector res = -instance.c;
+  for (HighsInt r = 0; r <res.dim; r++) {
+    for (HighsInt j = 0; j < r; j++) {
+      res.value[r] -= res.value[j] * L[j * instance.num_var + r];
+    }
+    res.value[r] /= L[r * instance.num_var + r];
+  }
+
+  for (HighsInt i = res.dim - 1; i >= 0; i--) {
+    double sum = 0.0;
+    for (HighsInt j = res.dim - 1; j > i; j--) {
+      sum += res.value[j] * L[i * instance.num_var + j];
+    }
+    res.value[i] = (res.value[i] - sum) / L[i * instance.num_var + i];
+  }
+
+  // project solution to bounds and collect active bounds
+  Vector x0(instance.num_var);
+  Vector ra(instance.num_con);
+  std::vector<HighsInt> initialactive;
+  std::vector<HighsInt> initialinactive;
+  std::vector<BasisStatus> atlower;
+
+  for (int i=0; i<instance.num_var; i++) {
+    if (res.value[i] > 0.5/settings.hessianregularizationfactor 
+        && instance.var_up[i] == std::numeric_limits<double>::infinity()
+        && instance.c.value[i] < 0.0) {
+          modelstatus = QpModelStatus::UNBOUNDED;
+          return;
+    } else if (res.value[i] < 0.5/settings.hessianregularizationfactor 
+        && instance.var_lo[i] == std::numeric_limits<double>::infinity()
+        && instance.c.value[i] > 0.0) {
+          modelstatus = QpModelStatus::UNBOUNDED;
+          return;
+    } else if (res.value[i] <= instance.var_lo[i]) {
+        res.value[i] = instance.var_lo[i];
+        initialactive.push_back(i + instance.num_con);
+        atlower.push_back(BasisStatus::ActiveAtLower);
+    } else if (res.value[i] >= instance.var_up[i]) {
+        res.value[i] = instance.var_up[i];
+        initialactive.push_back(i + instance.num_con);
+        atlower.push_back(BasisStatus::ActiveAtUpper);
+    } else {
+        initialinactive.push_back(i + instance.num_con);
+    }
+    if (fabs(res.value[i]) > 10E-5) {
+      x0.value[i] = res.value[i];
+      x0.index[x0.num_nz++] = i;
+    }
+  }
+
+  // if no bounds are active, solution lies in the interior -> optimal
+  if (initialactive.size() == 0) {
+    modelstatus = QpModelStatus::OPTIMAL;
+  }
+
+  assert((HighsInt)(initialactive.size() + initialinactive.size()) ==
+         instance.num_var);
+
+  result.status = atlower;
+  result.active = initialactive;
+  result.inactive = initialinactive;
+  result.primal = x0;
+  result.rowact = ra;
+}
+
+#endif

src/qpsolver/gradient.hpp

@@ -11,17 +11,17 @@ class Gradient {
   bool uptodate;
   HighsInt numupdates = 0;
 
+ public:
+  Gradient(Runtime& rt)
+      : runtime(rt), gradient(Vector(rt.instance.num_var)), uptodate(false) {}
+
   void recompute() {
     runtime.instance.Q.vec_mat(runtime.primal, gradient);
     gradient += runtime.instance.c;
     uptodate = true;
     numupdates = 0;
   }
 
- public:
-  Gradient(Runtime& rt)
-      : runtime(rt), gradient(Vector(rt.instance.num_var)), uptodate(false) {}
-
   Vector& getGradient() {
     if (!uptodate ||
         numupdates >= runtime.settings.gradientrecomputefrequency) {

src/qpsolver/pricing.hpp

@@ -8,6 +8,7 @@ class Pricing {
   virtual HighsInt price(const Vector& x, const Vector& gradient) = 0;
   virtual void update_weights(const Vector& aq, const Vector& ep, HighsInt p,
                               HighsInt q) = 0;
+  virtual void recompute() = 0;
   virtual ~Pricing() {}
 };