From 17eab7d77ea21ed7533e5fb92e4263aa9c752876 Mon Sep 17 00:00:00 2001 From: "lperron@google.com" Date: Mon, 19 Sep 2011 08:58:03 +0000 Subject: [PATCH] estimate numerical problems in linear solver --- linear_solver/cbc_interface.cc | 4 + linear_solver/clp_interface.cc | 5 + linear_solver/glpk_interface.cc | 155 +++++++++++++++++++++++++++++++ linear_solver/linear_solver.cc | 6 +- linear_solver/linear_solver.h | 26 ++++++ linear_solver/linear_solver.swig | 1 + linear_solver/scip_interface.cc | 4 + 7 files changed, 200 insertions(+), 1 deletion(-) diff --git a/linear_solver/cbc_interface.cc b/linear_solver/cbc_interface.cc index cd7868b3ce..f595dcdade 100644 --- a/linear_solver/cbc_interface.cc +++ b/linear_solver/cbc_interface.cc @@ -161,6 +161,10 @@ class CBCInterface : public MPSolverInterface { return reinterpret_cast(&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. diff --git a/linear_solver/clp_interface.cc b/linear_solver/clp_interface.cc index 7fa0dea031..df36195265 100644 --- a/linear_solver/clp_interface.cc +++ b/linear_solver/clp_interface.cc @@ -123,6 +123,11 @@ class CLPInterface : public MPSolverInterface { return reinterpret_cast(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(); diff --git a/linear_solver/glpk_interface.cc b/linear_solver/glpk_interface.cc index 812c997b9d..ef864f9f5e 100644 --- a/linear_solver/glpk_interface.cc +++ b/linear_solver/glpk_interface.cc @@ -12,6 +12,7 @@ // limitations under the License. // +#include #include #include "base/hash.h" #include @@ -170,6 +171,8 @@ class GLPKInterface : public MPSolverInterface { return reinterpret_cast(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 row_scaling_factor(new double[num_rows + 1]); + scoped_array 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 values(new double[num_rows + 1]); + scoped_array 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 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) { diff --git a/linear_solver/linear_solver.cc b/linear_solver/linear_solver.cc index b923d872b3..8a515d8c5a 100644 --- a/linear_solver/linear_solver.cc +++ b/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 diff --git a/linear_solver/linear_solver.h b/linear_solver/linear_solver.h index cb2af6ce29..38c173d873 100644 --- a/linear_solver/linear_solver.h +++ b/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_; diff --git a/linear_solver/linear_solver.swig b/linear_solver/linear_solver.swig index 0c06ef0b67..9b0602fa09 100644 --- a/linear_solver/linear_solver.swig +++ b/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; diff --git a/linear_solver/scip_interface.cc b/linear_solver/scip_interface.cc index a7ffcd3227..324654da85 100644 --- a/linear_solver/scip_interface.cc +++ b/linear_solver/scip_interface.cc @@ -127,6 +127,10 @@ class SCIPInterface : public MPSolverInterface { return reinterpret_cast(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);