estimate numerical problems in linear solver

linear_solver/cbc_interface.cc

@@ -161,6 +161,10 @@ class CBCInterface : public MPSolverInterface {
     return reinterpret_cast<void*>(&osi_);
   }
 
+  virtual double ComputeExactConditionNumber() const {
+    LOG(FATAL) << "Condition number only available for continuous problems";
+  }
+
  private:
   // Reset best objective bound to +/- infinity depending on the
   // optimization direction.

linear_solver/clp_interface.cc

@@ -123,6 +123,11 @@ class CLPInterface : public MPSolverInterface {
     return reinterpret_cast<void*>(clp_.get());
   }
 
+  virtual double ComputeExactConditionNumber() const {
+    // CLP does not provide the necessary API to access the inverse basis.
+    LOG(FATAL) << "ComputeExactConditionNumber is not implemented in CLP.";
+  }
+
  private:
   // Create dummy variable to be able to create empty constraints.
   void CreateDummyVariableForEmptyConstraints();

linear_solver/glpk_interface.cc

@@ -12,6 +12,7 @@
 // limitations under the License.
 //
 
+#include <math.h>
 #include <stddef.h>
 #include "base/hash.h"
 #include <limits>
@@ -170,6 +171,8 @@ class GLPKInterface : public MPSolverInterface {
     return reinterpret_cast<void*>(lp_);
   }
 
+  virtual double ComputeExactConditionNumber() const;
+
  private:
   // Configure the solver's parameters.
   void ConfigureGLPKParameters(const MPSolverParameters& param);
@@ -190,6 +193,20 @@ class GLPKInterface : public MPSolverInterface {
   // Transforms basis status from GLPK integer code to MPSolver::BasisStatus.
   MPSolver::BasisStatus TransformGLPKBasisStatus(int glpk_basis_status) const;
 
+  // Computes the L1-norm of the current scaled basis.
+  // The L1-norm |A| is defined as max_j sum_i |a_ij|
+  // This method is available only for continuous problems.
+  double ComputeScaledBasisL1Norm(
+    int num_rows, int num_cols,
+    double* row_scaling_factor, double* column_scaling_factor) const;
+
+  // Computes the L1-norm of the inverse of the current scaled
+  // basis.
+  // This method is available only for continuous problems.
+  double ComputeInverseScaledBasisL1Norm(
+    int num_rows, int num_cols,
+    double* row_scaling_factor, double* column_scaling_factor) const;
+
   glp_prob* lp_;
   bool mip_;
 
@@ -755,6 +772,144 @@ void GLPKInterface::CheckBestObjectiveBoundExists() const {
   }
 }
 
+double GLPKInterface::ComputeExactConditionNumber() const {
+  CHECK(IsContinuous()) <<
+      "Condition number only available for continuous problems";
+  CheckSolutionIsSynchronized();
+  // Simplex is the only LP algorithm supported in the wrapper for
+  // GLPK, so when a solution exists, a basis exists.
+  CheckSolutionExists();
+  const int num_rows = glp_get_num_rows(lp_);
+  const int num_cols = glp_get_num_cols(lp_);
+  // GLPK indexes everything starting from 1 instead of 0.
+  scoped_array<double> row_scaling_factor(new double[num_rows + 1]);
+  scoped_array<double> column_scaling_factor(new double[num_cols + 1]);
+  for (int row = 1; row <= num_rows; ++row) {
+    row_scaling_factor[row] = glp_get_rii(lp_, row);
+  }
+  for (int col = 1; col <= num_cols; ++col) {
+    column_scaling_factor[col] = glp_get_sjj(lp_, col);
+  }
+  return
+      ComputeInverseScaledBasisL1Norm(
+          num_rows, num_cols,
+          row_scaling_factor.get(), column_scaling_factor.get()) *
+      ComputeScaledBasisL1Norm(
+          num_rows, num_cols,
+          row_scaling_factor.get(), column_scaling_factor.get());
+}
+
+double GLPKInterface::ComputeScaledBasisL1Norm(
+    int num_rows, int num_cols,
+    double* row_scaling_factor, double* column_scaling_factor) const {
+  double norm = 0.0;
+  scoped_array<double> values(new double[num_rows + 1]);
+  scoped_array<int> indices(new int[num_rows + 1]);
+  for (int col = 1; col <= num_cols; ++col) {
+    const int glpk_basis_status = glp_get_col_stat(lp_, col);
+    // Take into account only basic columns.
+    if (glpk_basis_status == GLP_BS) {
+      // Compute L1-norm of column 'col': sum_row |a_row,col|.
+      const int num_nz = glp_get_mat_col(lp_, col, indices.get(), values.get());
+      double column_norm = 0.0;
+      for (int k = 1; k <= num_nz; k++) {
+        column_norm += fabs(values[k] * row_scaling_factor[indices[k]]);
+      }
+      column_norm *= fabs(column_scaling_factor[col]);
+      // Compute max_col column_norm
+      norm = std::max(norm, column_norm);
+    }
+  }
+  // Slack variables.
+  for (int row = 1; row <= num_rows; ++row) {
+    const int glpk_basis_status = glp_get_row_stat(lp_, row);
+    // Take into account only basic slack variables.
+    if (glpk_basis_status == GLP_BS) {
+      // Only one non-zero coefficient: +/- 1.0 in the corresponding
+      // row. The row has a scaling coefficient but the slack variable
+      // is never scaled on top of that.
+      const double column_norm = fabs(row_scaling_factor[row]);
+      // Compute max_col column_norm
+      norm = std::max(norm, column_norm);
+    }
+  }
+  return norm;
+}
+
+double GLPKInterface::ComputeInverseScaledBasisL1Norm(
+    int num_rows, int num_cols,
+    double* row_scaling_factor, double* column_scaling_factor) const {
+  // Compute the LU factorization if it doesn't exist yet.
+  if (!glp_bf_exists(lp_)) {
+    const int factorize_status = glp_factorize(lp_);
+    switch (factorize_status) {
+      case GLP_EBADB: {
+        LOG(FATAL) << "Not able to factorize: error GLP_EBADB.";
+        break;
+      }
+      case GLP_ESING: {
+        LOG(WARNING)
+            << "Not able to factorize: "
+            << "the basis matrix is singular within the working precision.";
+        return MPSolver::infinity();
+      }
+      case GLP_ECOND: {
+        LOG(WARNING)
+            << "Not able to factorize: the basis matrix is ill-conditioned.";
+        return MPSolver::infinity();
+      }
+      default:
+        break;
+    }
+  }
+  scoped_array<double> right_hand_side(new double[num_rows + 1]);
+  double norm = 0.0;
+  // Iteratively solve B x = e_k, where e_k is the kth unit vector.
+  // The result of this computation is the kth column of B^-1.
+  // glp_ftran works on original matrix. Scale input and result to
+  // obtain the norm of the kth column in the inverse scaled
+  // matrix. See glp_ftran documentation in glpapi12.c for how the
+  // scaling is done: inv(B'') = inv(SB) * inv(B) * inv(R) where:
+  // o B'' is the scaled basis
+  // o B is the original basis
+  // o R is the diagonal row scaling matrix
+  // o SB consists of the basic columns of the augmented column
+  // scaling matrix (auxiliary variables then structural variables):
+  // S~ = diag(inv(R) | S).
+  for (int k = 1; k <= num_rows; ++k) {
+    for (int row = 1; row <= num_rows; ++row) {
+      right_hand_side[row] = 0.0;
+    }
+    right_hand_side[k] = 1.0;
+    // Multiply input by inv(R).
+    for (int row = 1; row <= num_rows; ++row) {
+      right_hand_side[row] /= row_scaling_factor[row];
+    }
+    glp_ftran(lp_, right_hand_side.get());
+    // glp_ftran stores the result in the same vector where the right
+    // hand side was provided.
+    // Multiply result by inv(SB).
+    for (int row = 1; row <= num_rows; ++row) {
+      const int k = glp_get_bhead(lp_, row);
+      if (k <= num_rows) {
+        // Auxiliary variable.
+        right_hand_side[row] *= row_scaling_factor[k];
+      } else {
+        // Structural variable.
+        right_hand_side[row] /= column_scaling_factor[k - num_rows];
+      }
+    }
+    // Compute sum_row |vector_row|.
+    double column_norm = 0.0;
+    for (int row = 1; row <= num_rows; ++row) {
+      column_norm += fabs(right_hand_side[row]);
+    }
+    // Compute max_col column_norm
+    norm = std::max(norm, column_norm);
+  }
+  return norm;
+}
+
 // ------ Parameters ------
 
 void GLPKInterface::ConfigureGLPKParameters(const MPSolverParameters& param) {

linear_solver/linear_solver.cc

@@ -781,6 +781,10 @@ int64 MPSolver::nodes() const {
   return interface_->nodes();
 }
 
+double MPSolver::ComputeExactConditionNumber() const {
+  return interface_->ComputeExactConditionNumber();
+}
+
 // ---------- MPSolverInterface ----------
 
 const int MPSolverInterface::kDummyVariableIndex = 0;
@@ -887,7 +891,7 @@ void MPSolverInterface::InvalidateSolutionSynchronization() {
 void MPSolverInterface::SetCommonParameters(const MPSolverParameters& param) {
   SetPrimalTolerance(param.GetDoubleParam(
       MPSolverParameters::PRIMAL_TOLERANCE));
-  SetPrimalTolerance(param.GetDoubleParam(MPSolverParameters::DUAL_TOLERANCE));
+  SetDualTolerance(param.GetDoubleParam(MPSolverParameters::DUAL_TOLERANCE));
   SetPresolveMode(param.GetIntegerParam(MPSolverParameters::PRESOLVE));
   // TODO(user): In the future, we could distinguish between the
   // algorithm to solve the root LP and the algorithm to solve node

linear_solver/linear_solver.h

@@ -324,6 +324,28 @@ class MPSolver {
   // SCIP*).
   void* underlying_solver();
 
+  // Computes the exact condition number of the current scaled basis:
+  // L1norm(B) * L1norm(inverse(B)), where B is the scaled basis.
+  // This method requires that a basis exists: it should be called
+  // after Solve. It is only available for continuous problems. It is
+  // implemented for GLPK but not CLP because CLP does not provide the
+  // API for doing it.
+  // The condition number measures how well the constraint matrix is
+  // conditioned and can be used to predict whether numerical issues
+  // will arise during the solve: the model is declared infeasible
+  // whereas it is feasible (or vice-versa), the solution obtained is
+  // not optimal or violates some constraints, the resolution is slow
+  // because of repeated singularities.
+  // The rule of thumb to interpret the condition number kappa is:
+  // o kappa <= 1e7: virtually no chance of numerical issues
+  // o 1e7 < kappa <= 1e10: small chance of numerical issues
+  // o 1e10 < kappa <= 1e13: medium chance of numerical issues
+  // o kappa > 1e13: high chance of numerical issues
+  // The computation of the condition number depends on the quality of
+  // the LU decomposition, so it is not very accurate when the matrix
+  // is ill conditioned.
+  double ComputeExactConditionNumber() const;
+
   friend class GLPKInterface;
   friend class CLPInterface;
   friend class CBCInterface;
@@ -791,6 +813,10 @@ class MPSolverInterface {
   // Returns the underlying solver.
   virtual void* underlying_solver() = 0;
 
+  // Computes exact condition number. Only available for continuous
+  // problems and only implemented in GLPK.
+  virtual double ComputeExactConditionNumber() const = 0;
+
   friend class MPSolver;
  protected:
   MPSolver* const solver_;

linear_solver/linear_solver.swig

@@ -403,6 +403,7 @@ namespace operations_research {
 %rename (checkNameValidity) MPSolver::CheckNameValidity;
 %rename (clear) MPSolver::Clear;
 %rename (clearObjective) MPSolver::ClearObjective;
+%rename (computeExactConditionNumber) MPSolver::ComputeExactConditionNumber;
 %rename (init) MPSolver::Init;
 %rename (isLpFormat) MPSolver::IsLPFormat;
 %rename (loadModel) MPSolver::LoadModel;

linear_solver/scip_interface.cc

@@ -127,6 +127,10 @@ class SCIPInterface : public MPSolverInterface {
     return reinterpret_cast<void*>(scip_);
   }
 
+  virtual double ComputeExactConditionNumber() const {
+    LOG(FATAL) << "Condition number only available for continuous problems";
+  }
+
  private:
   // Set all parameters in the underlying solver.
   virtual void SetParameters(const MPSolverParameters& param);