Google OR-Tools v9.14
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
set_cover_cft.h
Go to the documentation of this file.
1// Copyright 2025 Francesco Cavaliere
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
14#ifndef OR_TOOLS_ORTOOLS_SET_COVER_SET_COVER_CFT_H
15#define OR_TOOLS_ORTOOLS_SET_COVER_SET_COVER_CFT_H
16
17#include <absl/algorithm/container.h>
18#include <absl/base/internal/pretty_function.h>
19#include <absl/status/status.h>
20
21#include <limits>
22
26
28
29// Implementation of:
30// Caprara, Alberto, Matteo Fischetti, and Paolo Toth. 1999. “A Heuristic
31// Method for the Set Covering Problem.” Operations Research 47 (5): 730–43.
32// https://www.jstor.org/stable/223097
33//
34// Hereafter referred to as CFT.
35//
36// SUMMARY
37// The CFT algorithm is a heuristic approach to the set-covering problem. At its
38// core, it combines a primal greedy heuristic with dual information obtained
39// from the optimization of the Lagrangian relaxation of the problem.
40// (Note: columns/subsets and rows/elements will be used interchangeably.)
41//
42// STRUCTURE
43// The core of the algorithm is the 3Phase:
44// 1. Subgradient optimization of the Lagrangian relaxation.
45// 2. A primal greedy heuristic guided by the dual information.
46// 3. Fixing some of the "best" columns (in terms of reduced costs) into the
47// solution (diving).
48// + Repeat until an exit criterion is met.
49//
50// The paper also considers an optional outer loop, which invokes the 3Phase
51// process and then fixes some columns from the current best solution. This
52// introduces two levels of diving: the outer loop fixes "primal-good" columns
53// (based on the best solution), while the inner loop fixes "dual-good" columns
54// (based on reduced costs).
55// NOTE: The outer loop is not implemented in this version (yet - April 2025).
56//
57// Key characteristics of the algorithm:
58//
59// - The CFT algorithm is tailored for instances where the number of columns is
60// significantly larger than the number of rows.
61//
62// - To improve efficiency, a core model approach is used. This involves
63// selecting a small subset of columns based on their reduced costs, thereby
64// substantially reducing the problem size handled in the internal steps.
65//
66// - Due to the use of the core model and column fixing, the algorithm rarely
67// considers the entire problem. Instead, it operates on a small "window" of
68// the problem. Efficiently managing this small window is a central aspect of
69// any CFT implementation.
70//
71// - The core model scheme also enables an alternative implementation where the
72// algorithm starts with a small model and progressively adds columns through
73// a column-generation procedure. While this column generation procedure is
74// problem-dependent and cannot be implemented here, the architecture of this
75// implementation is designed to be extensible, allowing for such a procedure
76// to be added in the future.
77//
78
82
83// Small class to store the solution of a sub-model. It contains the cost and
84// the subset list.
85class Solution {
86 public:
87 Solution() = default;
88 Solution(const SubModel& model, const std::vector<SubsetIndex>& core_subsets);
90 double cost() const { return cost_; }
91 const std::vector<FullSubsetIndex>& subsets() const { return subsets_; }
92 void AddSubset(FullSubsetIndex subset, Cost cost) {
93 subsets_.push_back(subset);
94 cost_ += cost;
95 }
96 bool Empty() const { return subsets_.empty(); }
97 void Clear() {
98 cost_ = 0.0;
99 subsets_.clear();
100 }
101
102 private:
103 Cost cost_ = std::numeric_limits<Cost>::max();
104 std::vector<FullSubsetIndex> subsets_;
105};
106
107// In the narrow scope of the CFT subgradient, there are often divisions
108// between non-negative quantities (e.g., to compute a relative gap). In these
109// specific cases, the denominator should always be greater than the
110// numerator. This function checks that.
111inline Cost DivideIfGE0(Cost numerator, Cost denominator) {
112 DCHECK_GE(numerator, -1e-6);
113 if (numerator < 1e-6) {
114 return 0.0;
115 }
116 return numerator / denominator;
117}
118
119// Dual information related to a SubModel.
120// Stores multipliers, reduced costs, and the lower bound, and provides an
121// interface that keeps them alligned.
122class DualState {
123 public:
124 DualState() = default;
125 DualState(const DualState&) = default;
126 template <typename SubModelT>
127 explicit DualState(const SubModelT& model)
128 : lower_bound_(.0),
129 multipliers_(model.num_elements(), .0),
130 reduced_costs_(model.subset_costs().begin(),
131 model.subset_costs().end()) {}
132
133 Cost lower_bound() const { return lower_bound_; }
134 const ElementCostVector& multipliers() const { return multipliers_; }
135 const SubsetCostVector& reduced_costs() const { return reduced_costs_; }
137 // NOTE: This function contains one of the two O(nnz) subgradient steps
138 template <typename SubModelT, typename Op>
139 void DualUpdate(const SubModelT& model, Op multiplier_operator) {
140 multipliers_.resize(model.num_elements());
141 reduced_costs_.resize(model.num_subsets());
142 lower_bound_ = .0;
143 // Update multipliers
144 for (ElementIndex i : model.ElementRange()) {
145 multiplier_operator(i, multipliers_[i]);
146 lower_bound_ += multipliers_[i];
147 DCHECK(std::isfinite(multipliers_[i]));
148 DCHECK_GE(multipliers_[i], .0);
149 }
150 lower_bound_ += ComputeReducedCosts(model, multipliers_, reduced_costs_);
151 }
152
153 private:
154 // Single hot point to optimize once for the different use cases.
155 template <typename ModelT>
156 static Cost ComputeReducedCosts(const ModelT& model,
159 // Compute new reduced costs (O(nnz))
160 Cost negative_sum = .0;
161 for (SubsetIndex j : model.SubsetRange()) {
162 reduced_costs[j] = model.subset_costs()[j];
163 for (ElementIndex i : model.columns()[j]) {
164 reduced_costs[j] -= multipliers[i];
165 }
166 if (reduced_costs[j] < .0) {
167 negative_sum += reduced_costs[j];
168 }
169 }
170 return negative_sum;
171 }
172
173 Cost lower_bound_;
174 ElementCostVector multipliers_;
175 SubsetCostVector reduced_costs_;
176};
177
178// Utility aggregate to store and pass around both primal and dual states.
179struct PrimalDualState {
180 Solution solution;
185
187
188// Utilitiy aggregate used by the SubgradientOptimization procedure to
189// communicate pass the needed information to the SubgradientCBs interface.
190struct SubgradientContext {
191 const SubModel& model;
193
194 // Avoid copying unused reduced cost during subgradient
197
198 const Solution& best_solution;
201
202// Generic set of callbacks hooks used to specialized the behavior of the
203// subgradient optimization
204class SubgradientCBs {
205 public:
206 virtual bool ExitCondition(const SubgradientContext&) = 0;
207 virtual void RunHeuristic(const SubgradientContext&, Solution&) = 0;
209 ElementCostVector& delta_mults) = 0;
211 CoreModel& core_model, bool force = false) = 0;
212 virtual ~SubgradientCBs() = default;
213};
215// Subgradient optimization procedure. Optimizes the Lagrangian relaxation of
216// the Set-Covering problem until a termination criterion si met.
217void SubgradientOptimization(SubModel& core_model, SubgradientCBs& cbs,
218 PrimalDualState& best_state);
219
220// Subgradient callbacks implementation focused on improving the current best
221// dual bound.
222class BoundCBs : public SubgradientCBs {
223 public:
224 static constexpr Cost kTol = 1e-6;
225
226 BoundCBs(const SubModel& model);
227 Cost step_size() const { return step_size_; }
228 bool ExitCondition(const SubgradientContext& context) override;
230 ElementCostVector& delta_mults) override;
231 void RunHeuristic(const SubgradientContext& context,
232 Solution& solution) override {}
233 bool UpdateCoreModel(SubgradientContext context, CoreModel& core_model,
234 bool force = false) override;
236 private:
237 void MakeMinimalCoverageSubgradient(const SubgradientContext& context,
238 ElementCostVector& subgradient);
239
240 private:
241 Cost squared_norm_;
242 ElementCostVector direction_;
243 std::vector<SubsetIndex> lagrangian_solution_;
244
245 // Stopping condition
246 Cost prev_best_lb_;
247 BaseInt max_iter_countdown_;
248 BaseInt exit_test_countdown_;
249 BaseInt exit_test_period_;
250 BaseInt unfixed_run_extension_;
251
252 // Step size
253 void UpdateStepSize(SubgradientContext context);
254 Cost step_size_;
255 Cost last_min_lb_seen_;
256 Cost last_max_lb_seen_;
257 BaseInt step_size_update_countdown_;
258 BaseInt step_size_update_period_;
259};
260
264
265// Higher level function that return a function obtained by the dual multiplier
266// based greedy.
268 const SubModel& model, const DualState& dual_state,
269 Cost cost_cutoff = std::numeric_limits<BaseInt>::max());
270
271// Lower level greedy function which creates or completes a "solution" (seen as
272// set of subsets of the current SubModel) until a cutoff size or cost is
273// reached.
274// Note: since the cutoff might be reached, the returned solution might not be
275// feasible.
276Cost CoverGreedly(const SubModel& model, const DualState& dual_state,
277 Cost cost_cutoff, BaseInt size_cutoff,
278 std::vector<SubsetIndex>& sol_subsets);
279
283
284// Subgradient callbacks implementation focused on wandering near the optimal
285// multipliers and invoke the multipliers based greedy heuristic at each
286// iteration.
287class HeuristicCBs : public SubgradientCBs {
288 public:
289 HeuristicCBs() : step_size_(0.1), countdown_(250) {};
290 void set_step_size(Cost step_size) { step_size_ = step_size; }
291 bool ExitCondition(const SubgradientContext& context) override {
292 Cost upper_bound =
293 context.best_solution.cost() - context.model.fixed_cost();
294 Cost lower_bound = context.best_lower_bound;
295 return upper_bound - .999 < lower_bound || --countdown_ <= 0;
296 }
297 void RunHeuristic(const SubgradientContext& context,
298 Solution& solution) override;
300 ElementCostVector& delta_mults) override;
301 bool UpdateCoreModel(SubgradientContext context, CoreModel& core_model,
302 bool force = false) override {
303 return false;
304 }
305
306 private:
307 Cost step_size_;
308 BaseInt countdown_;
310
312 const Solution& init_solution = {});
313
317
318PrimalDualState RunCftHeuristic(SubModel& model,
319 const Solution& init_solution = {});
320
324
325// Coverage counter to decide the number of columns to keep in the core model.
326static constexpr BaseInt kMinCov = 5;
327
328// CoreModel extractor. Stores a pointer to the full model and specilized
329// UpdateCore in such a way to updated the SubModel (stored as base class) and
330// focus the search on a small windows of the full model.
331class FullToCoreModel : public SubModel {
332 using base = SubModel;
333 struct UpdateTrigger {
334 BaseInt countdown;
335 BaseInt period;
336 BaseInt max_period;
337 };
339 // This class handles the logic for selecting columns based on their reduced
340 // costs and the number of rows they cover. While this implementation is more
341 // complex than what would typically be required for static `SetCoverModel`s,
342 // it is designed to efficiently handle dynamically updated models where new
343 // columns are generated over time. In this case, recomputing the row view
344 // from scratch each time would introduce significant overhead. To avoid this,
345 // the column selection logic operates solely on the column view, without
346 // relying on the row view.
347 //
348 // NOTE: A cleaner alternative would involve modifying the `SetCoverModel`
349 // implementation to support incremental updates to the row view as new
350 // columns are added. This approach would reduce overhead while enabling a
351 // simpler and more efficient column selection process.
352 //
353 // NOTE: The row-view based approach is available at commit:
354 // a598cf83930629853f72b964ebcff01f7a9378e0
355 class ColumnSelector {
356 public:
357 const std::vector<FullSubsetIndex>& ComputeNewSelection(
358 FilterModelView full_model,
359 const std::vector<FullSubsetIndex>& forced_columns,
360 const SubsetCostVector& reduced_costs);
361
362 private:
363 bool SelectColumn(FilterModelView full_model, SubsetIndex j);
364 void SelecteMinRedCostColumns(FilterModelView full_model,
365 const SubsetCostVector& reduced_costs);
366 void SelectMinRedCostByRow(FilterModelView full_model,
367 const SubsetCostVector& reduced_costs);
368
369 private:
370 std::vector<SubsetIndex> candidates_;
371 std::vector<SubsetIndex>::const_iterator first_unselected_;
372 ElementToIntVector row_cover_counts_;
373 BaseInt rows_left_to_cover_;
374
375 std::vector<FullSubsetIndex> selection_;
376 SubsetBoolVector selected_;
377 };
378
379 public:
380 FullToCoreModel() = default;
381 FullToCoreModel(const Model* full_model);
382 Cost FixMoreColumns(const std::vector<SubsetIndex>& columns_to_fix) override;
383 void ResetColumnFixing(const std::vector<FullSubsetIndex>& columns_to_fix,
384 const DualState& state) override;
385 bool UpdateCore(Cost best_lower_bound,
386 const ElementCostVector& best_multipliers,
387 const Solution& best_solution, bool force) override;
388 void ResetPricingPeriod();
389 const DualState& best_dual_state() const { return best_dual_state_; }
391
392 protected:
393 void SizeUpdate();
394 bool IsTimeToUpdate(Cost best_lower_bound, bool force);
395
396 decltype(auto) IsFocusCol(FullSubsetIndex j) {
397 return is_focus_col_[static_cast<SubsetIndex>(j)];
398 }
399 decltype(auto) IsFocusRow(FullElementIndex i) {
400 return is_focus_row_[static_cast<ElementIndex>(i)];
402 void UpdatePricingPeriod(const DualState& full_dual_state,
403 Cost core_lower_bound, Cost core_upper_bound);
404 Cost UpdateMultipliers(const ElementCostVector& core_multipliers);
405 void ComputeAndSetFocus(Cost best_lower_bound, const Solution& best_solution);
406
407 // Views are not composable (for now), so we can either access the full_model
408 // with the strongly typed view or with the filtered view.
409
410 // Access the full model filtered by the current columns fixed.
412 return FilterModelView(full_model_, &is_focus_col_, &is_focus_row_,
413 full_model_->num_subsets(),
414 full_model_->num_elements());
415 }
416
417 // Access the full model with the strongly typed view.
419 return StrongModelView(full_model_);
420 }
421
422 std::vector<FullSubsetIndex> SelectNewCoreColumns(
423 const std::vector<FullSubsetIndex>& forced_columns = {});
424
425 private:
426 const Model* full_model_;
427
428 // Note: The `is_focus_col_` vector duplicates information already present in
429 // `SubModelView::cols_sizes_`. However, it does not overlap with any data
430 // stored in `CoreModel`. Since `CoreModel` is expected to be the primary use
431 // case, this vector is explicitly maintained here to ensure compatibility.
432 SubsetBoolVector is_focus_col_;
433
434 // Note: The `is_focus_row_` vector is functionally redundant with either
435 // `CoreModel::full2core_row_map_` or `SubModelView::rows_sizes_`. These
436 // existing structures could be used to create the filtered view of the full
437 // model. However, doing so would require generalizing the current view
438 // system to work with generic functors instead of vectors of integral types.
439 // Since the number of elements is assumed to be not prohibitive, a simpler
440 // implementation that avoids this memory optimization was preferred.
441 ElementBoolVector is_focus_row_;
442
443 BaseInt selection_coefficient_ = kMinCov;
444 Cost prev_best_lower_bound_;
445 DualState full_dual_state_;
446 DualState best_dual_state_;
447
448 BaseInt update_countdown_;
449 BaseInt update_period_;
450 BaseInt update_max_period_;
451
452 ColumnSelector col_selector_; // Here to avoid reallocations
453};
454
455} // namespace operations_research::scp
456
457#endif /* OR_TOOLS_ORTOOLS_SET_COVER_SET_COVER_CFT_H */
void ComputeMultipliersDelta(const SubgradientContext &context, ElementCostVector &delta_mults) override
void RunHeuristic(const SubgradientContext &context, Solution &solution) override
bool ExitCondition(const SubgradientContext &context) override
bool UpdateCoreModel(SubgradientContext context, CoreModel &core_model, bool force=false) override
void DualUpdate(const SubModelT &model, Op multiplier_operator)
const SubsetCostVector & reduced_costs() const
const ElementCostVector & multipliers() const
DualState(const DualState &)=default
void ComputeAndSetFocus(Cost best_lower_bound, const Solution &best_solution)
std::vector< FullSubsetIndex > SelectNewCoreColumns(const std::vector< FullSubsetIndex > &forced_columns={})
bool UpdateCore(Cost best_lower_bound, const ElementCostVector &best_multipliers, const Solution &best_solution, bool force) override
const DualState & best_dual_state() const
FilterModelView FixingFullModelView() const
Access the full model filtered by the current columns fixed.
void UpdatePricingPeriod(const DualState &full_dual_state, Cost core_lower_bound, Cost core_upper_bound)
void ResetColumnFixing(const std::vector< FullSubsetIndex > &columns_to_fix, const DualState &state) override
Cost FixMoreColumns(const std::vector< SubsetIndex > &columns_to_fix) override
Cost UpdateMultipliers(const ElementCostVector &core_multipliers)
decltype(auto) IsFocusCol(FullSubsetIndex j)
StrongModelView StrongTypedFullModelView() const
Access the full model with the strongly typed view.
decltype(auto) IsFocusRow(FullElementIndex i)
bool IsTimeToUpdate(Cost best_lower_bound, bool force)
bool ExitCondition(const SubgradientContext &context) override
void RunHeuristic(const SubgradientContext &context, Solution &solution) override
void ComputeMultipliersDelta(const SubgradientContext &context, ElementCostVector &delta_mults) override
bool UpdateCoreModel(SubgradientContext context, CoreModel &core_model, bool force=false) override
void AddSubset(FullSubsetIndex subset, Cost cost)
const std::vector< FullSubsetIndex > & subsets() const
virtual void ComputeMultipliersDelta(const SubgradientContext &, ElementCostVector &delta_mults)=0
virtual void RunHeuristic(const SubgradientContext &, Solution &)=0
virtual bool ExitCondition(const SubgradientContext &)=0
virtual bool UpdateCoreModel(SubgradientContext context, CoreModel &core_model, bool force=false)=0
void resize(size_type new_size)
PrimalDualState RunThreePhase(SubModel &model, const Solution &init_solution)
PrimalDualState RunCftHeuristic(SubModel &model, const Solution &init_solution)
Cost CoverGreedly(const SubModel &model, const DualState &dual_state, Cost cost_cutoff, BaseInt stop_size, std::vector< SubsetIndex > &sol_subsets)
Cost DivideIfGE0(Cost numerator, Cost denominator)
static constexpr BaseInt kMinCov
Coverage counter to decide the number of columns to keep in the core model.
Solution RunMultiplierBasedGreedy(const SubModel &model, const DualState &dual_state, Cost cost_cutoff)
void SubgradientOptimization(SubModel &model, SubgradientCBs &cbs, PrimalDualState &best_state)
util_intops::StrongVector< SubsetIndex, Cost > SubsetCostVector
Definition base_types.h:58
Select next search node to expand Select next item_i to add this new search node to the search Generate a new search node where item_i is not in the knapsack Check validity of this new partial solution(using propagators) - If valid
util_intops::StrongVector< ElementIndex, Cost > ElementCostVector
Definition base_types.h:59
ClosedInterval::Iterator end(ClosedInterval interval)
util_intops::StrongVector< ElementIndex, BaseInt > ElementToIntVector
Definition base_types.h:64
util_intops::StrongVector< SubsetIndex, bool > SubsetBoolVector
Definition base_types.h:71
util_intops::StrongVector< ElementIndex, bool > ElementBoolVector
Definition base_types.h:72
double Cost
Basic non-strict type for cost. The speed penalty for using double is ~2%.
Definition base_types.h:32
ClosedInterval::Iterator begin(ClosedInterval interval)
Utility aggregate to store and pass around both primal and dual states.
const Cost & best_lower_bound
Avoid copying unused reduced cost during subgradient.