diff --git a/ortools/glop/basis_representation.cc b/ortools/glop/basis_representation.cc index ae152afd6d..a4f5553b99 100644 --- a/ortools/glop/basis_representation.cc +++ b/ortools/glop/basis_representation.cc @@ -210,6 +210,15 @@ Status BasisFactorization::Initialize() { return ComputeFactorization(); } +RowToColMapping BasisFactorization::ComputeInitialBasis( + const std::vector& 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() { diff --git a/ortools/glop/basis_representation.h b/ortools/glop/basis_representation.h index 0361c73ee7..34411f249f 100644 --- a/ortools/glop/basis_representation.h +++ b/ortools/glop/basis_representation.h @@ -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& candidates); + // Return the number of rows in the basis. RowIndex GetNumberOfRows() const { return compact_matrix_.num_rows(); } diff --git a/ortools/glop/lu_factorization.cc b/ortools/glop/lu_factorization.cc index d9aaa6ea41..fc937bc96a 100644 --- a/ortools/glop/lu_factorization.cc +++ b/ortools/glop/lu_factorization.cc @@ -65,6 +65,33 @@ Status LuFactorization::ComputeFactorization( return Status::OK(); } +RowToColMapping LuFactorization::ComputeInitialBasis( + const CompactSparseMatrix& matrix, + const std::vector& 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(); } diff --git a/ortools/glop/lu_factorization.h b/ortools/glop/lu_factorization.h index 58ef1efa72..3486dfa4f4 100644 --- a/ortools/glop/lu_factorization.h +++ b/ortools/glop/lu_factorization.h @@ -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& candidates); + // Returns the column permutation used by the LU factorization. const ColumnPermutation& GetColumnPermutation() const { return col_perm_; } diff --git a/ortools/glop/markowitz.cc b/ortools/glop/markowitz.cc index 28aee31139..c00cda3157 100644 --- a/ortools/glop/markowitz.cc +++ b/ortools/glop/markowitz.cc @@ -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(*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(*index) / - num_cols.value()); + stats_.basis_residual_singleton_column_ratio.Add( + static_cast(*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(); diff --git a/ortools/glop/revised_simplex.cc b/ortools/glop/revised_simplex.cc index e33dbfc82a..a9a9aa1184 100644 --- a/ortools/glop/revised_simplex.cc +++ b/ortools/glop/revised_simplex.cc @@ -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 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); } diff --git a/ortools/glop/variables_info.cc b/ortools/glop/variables_info.cc index 0a64224c5d..92d187f02b 100644 --- a/ortools/glop/variables_info.cc +++ b/ortools/glop/variables_info.cc @@ -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(); diff --git a/ortools/glop/variables_info.h b/ortools/glop/variables_info.h index 413f4b1824..52575dc544 100644 --- a/ortools/glop/variables_info.h +++ b/ortools/glop/variables_info.h @@ -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(); diff --git a/ortools/lp_data/sparse.h b/ortools/lp_data/sparse.h index ad53cc614b..6b9d491ffe 100644 --- a/ortools/lp_data/sparse.h +++ b/ortools/lp_data/sparse.h @@ -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* 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 columns_; }; // Specialization of a CompactSparseMatrix used for triangular matrices.