always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

make glop more robust to incorrect/incomplete basis

Browse Source
This commit is contained in:
Laurent Perron
2021-08-20 13:37:19 +02:00
parent 692867f460
commit 23c1f491de
9 changed files with 140 additions and 50 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
ortools/glop/basis_representation.cc
View File

@@ -210,6 +210,15 @@ Status BasisFactorization::Initialize() {
return ComputeFactorization();
}
RowToColMapping BasisFactorization::ComputeInitialBasis(
const std::vector<ColIndex>& candidates) {
const RowToColMapping basis =
lu_factorization_.ComputeInitialBasis(compact_matrix_, candidates);
deterministic_time_ +=
lu_factorization_.DeterministicTimeOfLastFactorization();
return basis;
}
bool BasisFactorization::IsRefactorized() const { return num_updates_ == 0; }
Status BasisFactorization::Refactorize() {

12
ortools/glop/basis_representation.h
View File

@@ -188,6 +188,18 @@ class BasisFactorization {
// could not be factorized.
ABSL_MUST_USE_RESULT Status Initialize();
// This mainly forward the call to LuFactorization::ComputeInitialBasis().
//
// Note that once this is called, one would need to call Initialize() to
// actually create the factorization. The only side effect of this is to
// update the deterministic time.
//
// TODO(user): This "double" factorization is a bit inefficient, and we should
// probably Initialize() right away the factorization with the new basis, but
// more code is needed for that. It is also not that easy also because we want
// to permute all the added slack first.
RowToColMapping ComputeInitialBasis(const std::vector<ColIndex>& candidates);
// Return the number of rows in the basis.
RowIndex GetNumberOfRows() const { return compact_matrix_.num_rows(); }

27
ortools/glop/lu_factorization.cc
View File

@@ -65,6 +65,33 @@ Status LuFactorization::ComputeFactorization(
return Status::OK();
}
RowToColMapping LuFactorization::ComputeInitialBasis(
const CompactSparseMatrix& matrix,
const std::vector<ColIndex>& candidates) {
CompactSparseMatrixView view(&matrix, &candidates);
(void)markowitz_.ComputeRowAndColumnPermutation(view, &row_perm_, &col_perm_);
// Starts by the missing slacks.
RowToColMapping basis;
for (RowIndex row(0); row < matrix.num_rows(); ++row) {
if (row_perm_[row] == kInvalidRow) {
// Add the slack for this row.
basis.push_back(matrix.num_cols() +
RowToColIndex(row - matrix.num_rows()));
}
}
// Then add the used candidate columns.
CHECK_EQ(col_perm_.size(), candidates.size());
for (int i = 0; i < col_perm_.size(); ++i) {
if (col_perm_[ColIndex(i)] != kInvalidCol) {
basis.push_back(candidates[i]);
}
}
return basis;
}
double LuFactorization::DeterministicTimeOfLastFactorization() const {
return markowitz_.DeterministicTimeOfLastFactorization();
}

7
ortools/glop/lu_factorization.h
View File

@@ -55,6 +55,13 @@ class LuFactorization {
ABSL_MUST_USE_RESULT Status
ComputeFactorization(const CompactSparseMatrixView& compact_matrix);
// Given a set of columns, find a maximum linearly independent subset that can
// be factorized in a stable way, and complete it into a square matrix using
// slack columns. The initial set can have less, more or the same number of
// columns as the number of rows.
RowToColMapping ComputeInitialBasis(const CompactSparseMatrix& matrix,
const std::vector<ColIndex>& candidates);
// Returns the column permutation used by the LU factorization.
const ColumnPermutation& GetColumnPermutation() const { return col_perm_; }

22
ortools/glop/markowitz.cc
View File

@@ -140,9 +140,9 @@ Status Markowitz::ComputeRowAndColumnPermutation(
num_fp_operations_ += 10 * upper_.num_entries().value();
stats_.pivots_without_fill_in_ratio.Add(
1.0 * stats_num_pivots_without_fill_in / end_index);
1.0 * stats_num_pivots_without_fill_in / num_rows.value());
stats_.degree_two_pivot_columns.Add(1.0 * stats_degree_two_pivot_columns /
end_index);
num_rows.value());
return Status::OK();
}
@@ -218,7 +218,7 @@ void Markowitz::ExtractSingletonColumns(
}
}
stats_.basis_singleton_column_ratio.Add(static_cast<double>(*index) /
num_cols.value());
basis_matrix.num_rows().value());
}
bool Markowitz::IsResidualSingletonColumn(const ColumnView& column,
@@ -250,8 +250,8 @@ void Markowitz::ExtractResidualSingletonColumns(
upper_.AddTriangularColumn(column, row);
++(*index);
}
stats_.basis_residual_singleton_column_ratio.Add(static_cast<double>(*index) /
num_cols.value());
stats_.basis_residual_singleton_column_ratio.Add(
static_cast<double>(*index) / basis_matrix.num_rows().value());
}
const SparseColumn& Markowitz::ComputeColumn(const RowPermutation& row_perm,
@@ -817,12 +817,15 @@ void ColumnPriorityQueue::Reset(int max_degree, ColIndex num_cols) {
col_degree_.assign(num_cols, 0);
col_index_.assign(num_cols, -1);
col_by_degree_.resize(max_degree + 1);
min_degree_ = num_cols.value();
min_degree_ = max_degree + 1;
}
void ColumnPriorityQueue::PushOrAdjust(ColIndex col, int32_t degree) {
DCHECK_GE(degree, 0);
DCHECK_LT(degree, col_by_degree_.size());
DCHECK_GE(col, 0);
DCHECK_LT(col, col_degree_.size());
const int32_t old_degree = col_degree_[col];
if (degree != old_degree) {
const int32_t old_index = col_index_[col];
@@ -844,9 +847,12 @@ void ColumnPriorityQueue::PushOrAdjust(ColIndex col, int32_t degree) {
}
ColIndex ColumnPriorityQueue::Pop() {
while (col_by_degree_[min_degree_].empty()) {
min_degree_++;
DCHECK_GE(min_degree_, 0);
DCHECK_LE(min_degree_, col_by_degree_.size());
while (true) {
if (min_degree_ == col_by_degree_.size()) return kInvalidCol;
if (!col_by_degree_[min_degree_].empty()) break;
min_degree_++;
}
const ColIndex col = col_by_degree_[min_degree_].back();
col_by_degree_[min_degree_].pop_back();

61
ortools/glop/revised_simplex.cc
View File

@@ -68,6 +68,7 @@ class Cleanup {
DCHECK_GT(num_cols_, col); \
}
// TODO(user): Remove this function.
#define DCHECK_ROW_BOUNDS(row) \
{ \
DCHECK_LE(0, row); \
@@ -1353,32 +1354,51 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) {
// reuse the variable statuses.
const bool log_info = parameters_.log_search_progress() || VLOG_IS_ON(1);
if (solve_from_scratch && !solution_state_.IsEmpty()) {
// If an external basis has been provided or if the matrix changed, we need
// to perform more work, e.g., factorize the proposed basis and validate it.
variables_info_.InitializeFromBasisState(first_slack_col_, ColIndex(0),
solution_state_);
basis_.assign(num_rows_, kInvalidCol);
RowIndex row(0);
for (ColIndex col : variables_info_.GetIsBasicBitRow()) {
basis_[row] = col;
++row;
}
basis_factorization_.Clear();
reduced_costs_.ClearAndRemoveCostShifts();
primal_edge_norms_.Clear();
dual_edge_norms_.Clear();
dual_pricing_vector_.clear();
// TODO(user): If the basis is incomplete, we could complete it with
// better slack variables than is done by InitializeFirstBasis() by
// using a partial LU decomposition (see markowitz.h).
if (InitializeFirstBasis(basis_).ok()) {
solve_from_scratch = false;
} else {
if (log_info) {
LOG(INFO) << "RevisedSimplex is not using the warm start "
"basis because it is not factorizable.";
// If an external basis has been provided or if the matrix changed, we need
// to perform more work, e.g., factorize the proposed basis and validate it.
variables_info_.InitializeFromBasisState(first_slack_col_, ColIndex(0),
solution_state_);
// Use the set of basic columns as a "hint" to construct the first basis.
std::vector<ColIndex> candidates;
for (const ColIndex col : variables_info_.GetIsBasicBitRow()) {
candidates.push_back(col);
}
// Optimization: Try to factorize it right away if we have the correct
// number of element. Ideally the other path below would no require a
// "double" factorization effort, so this would not be needed.
if (candidates.size() == num_rows_) {
basis_.clear();
for (const ColIndex col : candidates) {
basis_.push_back(col);
}
// TODO(user): Depending on the error here, there is no point doing extra
// work below. This is the case when we fail because of a bad initial
// condition number for instance.
if (InitializeFirstBasis(basis_).ok()) {
solve_from_scratch = false;
}
}
if (solve_from_scratch) {
basis_ = basis_factorization_.ComputeInitialBasis(candidates);
variables_info_.CorrectBasicStatus(basis_);
if (InitializeFirstBasis(basis_).ok()) {
solve_from_scratch = false;
} else {
if (log_info) {
LOG(INFO) << "RevisedSimplex is not using the warm start "
"basis because it is not factorizable.";
}
}
}
}
@@ -3080,7 +3100,6 @@ Status RevisedSimplex::DualMinimize(bool feasibility_phase,
}
ColIndex RevisedSimplex::SlackColIndex(RowIndex row) const {
// TODO(user): Remove this function.
DCHECK_ROW_BOUNDS(row);
return first_slack_col_ + RowToColIndex(row);
}

36
ortools/glop/variables_info.cc
View File

@@ -107,9 +107,7 @@ void VariablesInfo::InitializeFromBasisState(ColIndex first_slack_col,
const BasisState& state) {
ResetStatusInfo();
RowIndex num_basic_variables(0);
const ColIndex num_cols = lower_bounds_.size();
const RowIndex num_rows = ColToRowIndex(num_cols - first_slack_col);
DCHECK_LE(num_new_cols, first_slack_col);
const ColIndex first_new_col(first_slack_col - num_new_cols);
@@ -131,21 +129,12 @@ void VariablesInfo::InitializeFromBasisState(ColIndex first_slack_col,
// Remove incompatibilities between the warm status and the current state.
switch (status) {
case VariableStatus::BASIC:
// Do not allow more than num_rows VariableStatus::BASIC variables.
if (num_basic_variables == num_rows) {
VLOG(1) << "Too many basic variables in the warm-start basis."
<< "Only keeping the first ones as VariableStatus::BASIC.";
UpdateToNonBasicStatus(col, DefaultVariableStatus(col));
} else {
++num_basic_variables;
// Because we just called ResetStatusInfo(), we optimize the call to
// UpdateToNonBasicStatus(col) here. In an incremental setting with
// almost no work per call, the update of all the DenseBitRow are
// visible.
variable_status_[col] = VariableStatus::BASIC;
is_basic_.Set(col, true);
}
// Because we just called ResetStatusInfo(), we optimize the call to
// UpdateToNonBasicStatus(col) here. In an incremental setting with
// almost no work per call, the update of all the DenseBitRow are
// visible.
variable_status_[col] = VariableStatus::BASIC;
is_basic_.Set(col, true);
break;
case VariableStatus::AT_LOWER_BOUND:
if (lower_bounds_[col] == upper_bounds_[col]) {
@@ -171,6 +160,19 @@ void VariablesInfo::InitializeFromBasisState(ColIndex first_slack_col,
}
}
void VariablesInfo::CorrectBasicStatus(const RowToColMapping& basis) {
const ColIndex num_cols = lower_bounds_.size();
is_basic_.ClearAndResize(num_cols);
for (const ColIndex col : basis) {
UpdateToBasicStatus(col);
}
for (ColIndex col(0); col < num_cols; ++col) {
if (!is_basic_[col] && variable_status_[col] == VariableStatus::BASIC) {
UpdateToNonBasicStatus(col, DefaultVariableStatus(col));
}
}
}
void VariablesInfo::InitializeToDefaultStatus() {
ResetStatusInfo();
const ColIndex num_cols = lower_bounds_.size();

4
ortools/glop/variables_info.h
View File

@@ -80,6 +80,10 @@ class VariablesInfo {
void InitializeFromBasisState(ColIndex first_slack, ColIndex num_new_cols,
const BasisState& state);
// Reset to the default status any column not listed in the basis and make
// sure all the one listed are marked as basic.
void CorrectBasicStatus(const RowToColMapping& basis);
// Sets all variables status to their lowest magnitude bounds. Note that there
// will be no basic variable after this is called.
void InitializeToDefaultStatus();

12
ortools/lp_data/sparse.h
View File

@@ -477,14 +477,18 @@ class CompactSparseMatrixView {
public:
CompactSparseMatrixView(const CompactSparseMatrix* compact_matrix,
const RowToColMapping* basis)
: compact_matrix_(*compact_matrix), basis_(*basis) {}
: compact_matrix_(*compact_matrix),
columns_(basis->data(), basis->size().value()) {}
CompactSparseMatrixView(const CompactSparseMatrix* compact_matrix,
const std::vector<ColIndex>* columns)
: compact_matrix_(*compact_matrix), columns_(*columns) {}
// Same behavior as the SparseMatrix functions above.
bool IsEmpty() const { return compact_matrix_.IsEmpty(); }
RowIndex num_rows() const { return compact_matrix_.num_rows(); }
ColIndex num_cols() const { return RowToColIndex(basis_.size()); }
ColIndex num_cols() const { return ColIndex(columns_.size()); }
const ColumnView column(ColIndex col) const {
return compact_matrix_.column(basis_[ColToRowIndex(col)]);
return compact_matrix_.column(columns_[col.value()]);
}
EntryIndex num_entries() const;
Fractional ComputeOneNorm() const;
@@ -494,7 +498,7 @@ class CompactSparseMatrixView {
// We require that the underlying CompactSparseMatrix and RowToColMapping
// continue to own the (potentially large) data accessed via this view.
const CompactSparseMatrix& compact_matrix_;
const RowToColMapping& basis_;
const absl::Span<const ColIndex> columns_;
};
// Specialization of a CompactSparseMatrix used for triangular matrices.