OR-Tools  8.2
initial_basis.cc
Go to the documentation of this file.
1// Copyright 2010-2018 Google LLC
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
15
16#include <queue>
17
20
21namespace operations_research {
22namespace glop {
23
25 const DenseRow& objective,
26 const DenseRow& lower_bound,
27 const DenseRow& upper_bound,
28 const VariableTypeRow& variable_type)
29 : max_scaled_abs_cost_(0.0),
30 bixby_column_comparator_(*this),
31 triangular_column_comparator_(*this),
32 compact_matrix_(compact_matrix),
33 objective_(objective),
34 lower_bound_(lower_bound),
35 upper_bound_(upper_bound),
36 variable_type_(variable_type) {}
37
38void InitialBasis::CompleteBixbyBasis(ColIndex num_cols,
39 RowToColMapping* basis) {
40 // Initialize can_be_replaced ('I' in Bixby's paper) and has_zero_coefficient
41 // ('r' in Bixby's paper).
42 const RowIndex num_rows = compact_matrix_.num_rows();
43 DenseBooleanColumn can_be_replaced(num_rows, false);
44 DenseBooleanColumn has_zero_coefficient(num_rows, false);
45 DCHECK_EQ(num_rows, basis->size());
46 basis->resize(num_rows, kInvalidCol);
47 for (RowIndex row(0); row < num_rows; ++row) {
48 if ((*basis)[row] == kInvalidCol) {
49 can_be_replaced[row] = true;
50 has_zero_coefficient[row] = true;
51 }
52 }
53
54 // This is 'v' in Bixby's paper.
55 DenseColumn scaled_diagonal_abs(compact_matrix_.num_rows(), kInfinity);
56
57 // Compute a list of candidate indices and sort them using the heuristic
58 // described in Bixby's paper.
59 std::vector<ColIndex> candidates;
60 ComputeCandidates(num_cols, &candidates);
61
62 // Loop over the candidate columns, and add them to the basis if the
63 // heuristics are satisfied.
64 for (int i = 0; i < candidates.size(); ++i) {
65 bool enter_basis = false;
66 const ColIndex candidate_col_index = candidates[i];
67 const auto& candidate_col = compact_matrix_.column(candidate_col_index);
68
69 // Bixby's heuristic only works with scaled columns. This should be the
70 // case by default since we only use this when the matrix is scaled, but
71 // it is not the case for our tests... The overhead for computing the
72 // infinity norm for each column should be minimal.
73 if (InfinityNorm(candidate_col) != 1.0) continue;
74
75 RowIndex candidate_row;
76 Fractional candidate_coeff = RestrictedInfinityNorm(
77 candidate_col, has_zero_coefficient, &candidate_row);
78 const Fractional kBixbyHighThreshold = 0.99;
79 if (candidate_coeff > kBixbyHighThreshold) {
80 enter_basis = true;
81 } else if (IsDominated(candidate_col, scaled_diagonal_abs)) {
82 candidate_coeff = RestrictedInfinityNorm(candidate_col, can_be_replaced,
83 &candidate_row);
84 if (candidate_coeff != 0.0) {
85 enter_basis = true;
86 }
87 }
88
89 if (enter_basis) {
90 can_be_replaced[candidate_row] = false;
91 SetSupportToFalse(candidate_col, &has_zero_coefficient);
92 const Fractional kBixbyLowThreshold = 0.01;
93 scaled_diagonal_abs[candidate_row] =
94 kBixbyLowThreshold * std::abs(candidate_coeff);
95 (*basis)[candidate_row] = candidate_col_index;
96 }
97 }
98}
99
101 RowToColMapping* basis) {
102 return GetMarosBasis<false>(num_cols, basis);
103}
104
105void InitialBasis::GetDualMarosBasis(ColIndex num_cols,
106 RowToColMapping* basis) {
107 return GetMarosBasis<true>(num_cols, basis);
108}
109
111 RowToColMapping* basis) {
112 return CompleteTriangularBasis<false>(num_cols, basis);
113}
114
116 RowToColMapping* basis) {
117 return CompleteTriangularBasis<true>(num_cols, basis);
118}
119
120template <bool only_allow_zero_cost_column>
121void InitialBasis::CompleteTriangularBasis(ColIndex num_cols,
122 RowToColMapping* basis) {
123 // Initialize can_be_replaced.
124 const RowIndex num_rows = compact_matrix_.num_rows();
125 DenseBooleanColumn can_be_replaced(num_rows, false);
126 DCHECK_EQ(num_rows, basis->size());
127 basis->resize(num_rows, kInvalidCol);
128 for (RowIndex row(0); row < num_rows; ++row) {
129 if ((*basis)[row] == kInvalidCol) {
130 can_be_replaced[row] = true;
131 }
132 }
133
134 // Initialize the residual non-zero pattern for the rows that can be replaced.
135 MatrixNonZeroPattern residual_pattern;
136 residual_pattern.Reset(num_rows, num_cols);
137 for (ColIndex col(0); col < num_cols; ++col) {
138 if (only_allow_zero_cost_column && objective_[col] != 0.0) continue;
139 for (const SparseColumn::Entry e : compact_matrix_.column(col)) {
140 if (can_be_replaced[e.row()]) {
141 residual_pattern.AddEntry(e.row(), col);
142 }
143 }
144 }
145
146 // Initialize a priority queue of residual singleton columns.
147 // Also compute max_scaled_abs_cost_ for GetColumnPenalty().
148 std::vector<ColIndex> residual_singleton_column;
149 max_scaled_abs_cost_ = 0.0;
150 for (ColIndex col(0); col < num_cols; ++col) {
151 max_scaled_abs_cost_ =
152 std::max(max_scaled_abs_cost_, std::abs(objective_[col]));
153 if (residual_pattern.ColDegree(col) == 1) {
154 residual_singleton_column.push_back(col);
155 }
156 }
157 const Fractional kBixbyWeight = 1000.0;
158 max_scaled_abs_cost_ =
159 (max_scaled_abs_cost_ == 0.0) ? 1.0 : kBixbyWeight * max_scaled_abs_cost_;
160 std::priority_queue<ColIndex, std::vector<ColIndex>,
161 InitialBasis::TriangularColumnComparator>
162 queue(residual_singleton_column.begin(), residual_singleton_column.end(),
163 triangular_column_comparator_);
164
165 // Process the residual singleton columns by priority and add them to the
166 // basis if their "diagonal" coefficient is not too small.
167 while (!queue.empty()) {
168 const ColIndex candidate = queue.top();
169 queue.pop();
170 if (residual_pattern.ColDegree(candidate) != 1) continue;
171
172 // Find the position of the singleton and compute the infinity norm of
173 // the column (note that this is always 1.0 if the problem was scaled).
174 RowIndex row(kInvalidRow);
175 Fractional coeff = 0.0;
176 Fractional max_magnitude = 0.0;
177 for (const SparseColumn::Entry e : compact_matrix_.column(candidate)) {
178 max_magnitude = std::max(max_magnitude, std::abs(e.coefficient()));
179 if (can_be_replaced[e.row()]) {
180 row = e.row();
181 coeff = e.coefficient();
182 break;
183 }
184 }
185 const Fractional kStabilityThreshold = 0.01;
186 if (std::abs(coeff) < kStabilityThreshold * max_magnitude) continue;
188
189 // Use this candidate column in the basis.
190 (*basis)[row] = candidate;
191 can_be_replaced[row] = false;
192 residual_pattern.DeleteRowAndColumn(row, candidate);
193 for (const ColIndex col : residual_pattern.RowNonZero(row)) {
194 if (col == candidate) continue;
195 residual_pattern.DecreaseColDegree(col);
196 if (residual_pattern.ColDegree(col) == 1) {
197 queue.push(col);
198 }
199 }
200 }
201}
202
203int InitialBasis::GetMarosPriority(ColIndex col) const {
204 // Priority values for columns as defined in Maros's book.
205 switch (variable_type_[col]) {
207 return 3;
209 return 2;
211 return 2;
213 return 1;
215 return 0;
216 }
217}
218
219int InitialBasis::GetMarosPriority(RowIndex row) const {
220 // Priority values for rows are equal to
221 // 3 - row priority values as defined in Maros's book
222 ColIndex slack_index(RowToColIndex(row) + compact_matrix_.num_cols() -
223 RowToColIndex(compact_matrix_.num_rows()));
224
225 return GetMarosPriority(slack_index);
226}
227
228template <bool only_allow_zero_cost_column>
229void InitialBasis::GetMarosBasis(ColIndex num_cols, RowToColMapping* basis) {
230 VLOG(1) << "Starting Maros crash procedure.";
231
232 // Initialize basis to the all-slack basis.
233 const RowIndex num_rows = compact_matrix_.num_rows();
234 const ColIndex first_slack = num_cols - RowToColIndex(num_rows);
235 DCHECK_EQ(num_rows, basis->size());
236 basis->resize(num_rows);
237 for (RowIndex row(0); row < num_rows; row++) {
238 (*basis)[row] = first_slack + RowToColIndex(row);
239 }
240
241 // Initialize the set of available rows and columns.
242 DenseBooleanRow available(num_cols, true);
243 for (ColIndex col(0); col < first_slack; ++col) {
244 if (variable_type_[col] == VariableType::FIXED_VARIABLE ||
245 (only_allow_zero_cost_column && objective_[col] != 0.0)) {
246 available[col] = false;
247 }
248 }
249 for (ColIndex col = first_slack; col < num_cols; ++col) {
250 if (variable_type_[col] == VariableType::UNCONSTRAINED) {
251 available[col] = false;
252 }
253 }
254
255 // Initialize the residual non-zero pattern for the active part of the matrix.
256 MatrixNonZeroPattern residual_pattern;
257 residual_pattern.Reset(num_rows, num_cols);
258 for (ColIndex col(0); col < first_slack; ++col) {
259 for (const SparseColumn::Entry e : compact_matrix_.column(col)) {
260 if (available[RowToColIndex(e.row())] && available[col]) {
261 residual_pattern.AddEntry(e.row(), col);
262 }
263 }
264 }
265
266 // Go over residual pattern and mark rows as unavailable.
267 for (RowIndex row(0); row < num_rows; row++) {
268 if (residual_pattern.RowDegree(row) == 0) {
269 available[RowToColIndex(row) + first_slack] = false;
270 }
271 }
272
273 for (;;) {
274 // Make row selection by the Row Priority Function (RPF) from Maros's
275 // book.
276 int max_row_priority_function = std::numeric_limits<int>::min();
277 RowIndex max_rpf_row = kInvalidRow;
278 for (RowIndex row(0); row < num_rows; row++) {
279 if (available[RowToColIndex(row) + first_slack]) {
280 const int rpf =
281 10 * (3 - GetMarosPriority(row)) - residual_pattern.RowDegree(row);
282 if (rpf > max_row_priority_function) {
283 max_row_priority_function = rpf;
284 max_rpf_row = row;
285 }
286 }
287 }
288 if (max_rpf_row == kInvalidRow) break;
289
290 // Trace row for nonzero entries and pick one with best Column Priority
291 // Function (cpf).
292 const Fractional kStabilityThreshold = 1e-3;
293 ColIndex max_cpf_col(kInvalidCol);
294 int max_col_priority_function(std::numeric_limits<int>::min());
295 Fractional pivot_absolute_value = 0.0;
296 for (const ColIndex col : residual_pattern.RowNonZero(max_rpf_row)) {
297 if (!available[col]) continue;
298 const int cpf =
299 10 * GetMarosPriority(col) - residual_pattern.ColDegree(col);
300 if (cpf > max_col_priority_function) {
301 // Make sure that the pivotal entry is not too small in magnitude.
302 Fractional max_magnitude = 0;
303 pivot_absolute_value = 0.0;
304 const auto& column_values = compact_matrix_.column(col);
305 for (const SparseColumn::Entry e : column_values) {
306 const Fractional absolute_value = std::fabs(e.coefficient());
307 if (e.row() == max_rpf_row) pivot_absolute_value = absolute_value;
308 max_magnitude = std::max(max_magnitude, absolute_value);
309 }
310 if (pivot_absolute_value >= kStabilityThreshold * max_magnitude) {
311 max_col_priority_function = cpf;
312 max_cpf_col = col;
313 }
314 }
315 }
316
317 if (max_cpf_col == kInvalidCol) {
318 available[RowToColIndex(max_rpf_row) + first_slack] = false;
319 continue;
320 }
321
322 // Ensure that the row leaving the basis has a lower priority than the
323 // column entering the basis. If the best column is not good enough mark
324 // row as unavailable and choose another one.
325 const int row_priority = GetMarosPriority(max_rpf_row);
326 const int column_priority = GetMarosPriority(max_cpf_col);
327 if (row_priority >= column_priority) {
328 available[RowToColIndex(max_rpf_row) + first_slack] = false;
329 continue;
330 }
331
332 // Use this candidate column in the basis. Update residual pattern and row
333 // counts list.
334 (*basis)[max_rpf_row] = max_cpf_col;
335
336 VLOG(2) << "Slack variable " << max_rpf_row << " replaced by column "
337 << max_cpf_col
338 << ". Pivot coefficient magnitude: " << pivot_absolute_value << ".";
339
340 available[max_cpf_col] = false;
341 available[first_slack + RowToColIndex(max_rpf_row)] = false;
342
343 // Maintain the invariant that all the still available columns will have
344 // zeros on the rows we already replaced. This ensures the lower-triangular
345 // nature (after permutation) of the returned basis.
346 residual_pattern.DeleteRowAndColumn(max_rpf_row, max_cpf_col);
347 for (const ColIndex col : residual_pattern.RowNonZero(max_rpf_row)) {
348 available[col] = false;
349 }
350 }
351}
352
353void InitialBasis::ComputeCandidates(ColIndex num_cols,
354 std::vector<ColIndex>* candidates) {
355 candidates->clear();
356 max_scaled_abs_cost_ = 0.0;
357 for (ColIndex col(0); col < num_cols; ++col) {
358 if (variable_type_[col] != VariableType::FIXED_VARIABLE &&
359 compact_matrix_.column(col).num_entries() > 0) {
360 candidates->push_back(col);
361 max_scaled_abs_cost_ =
362 std::max(max_scaled_abs_cost_, std::abs(objective_[col]));
363 }
364 }
365 const Fractional kBixbyWeight = 1000.0;
366 max_scaled_abs_cost_ =
367 (max_scaled_abs_cost_ == 0.0) ? 1.0 : kBixbyWeight * max_scaled_abs_cost_;
368 std::sort(candidates->begin(), candidates->end(), bixby_column_comparator_);
369}
370
371int InitialBasis::GetColumnCategory(ColIndex col) const {
372 // Only the relative position of the returned number is important, so we use
373 // 2 for the category C2 in Bixby's paper and so on.
374 switch (variable_type_[col]) {
376 return 2;
378 return 3;
380 return 3;
382 return 4;
384 return 5;
385 }
386}
387
388Fractional InitialBasis::GetColumnPenalty(ColIndex col) const {
389 const VariableType type = variable_type_[col];
390 Fractional penalty = 0.0;
391 if (type == VariableType::LOWER_BOUNDED) {
392 penalty = lower_bound_[col];
393 }
394 if (type == VariableType::UPPER_BOUNDED) {
395 penalty = -upper_bound_[col];
396 }
398 penalty = lower_bound_[col] - upper_bound_[col];
399 }
400 return penalty + std::abs(objective_[col]) / max_scaled_abs_cost_;
401}
402
403bool InitialBasis::BixbyColumnComparator::operator()(ColIndex col_a,
404 ColIndex col_b) const {
405 if (col_a == col_b) return false;
406 const int category_a = initial_basis_.GetColumnCategory(col_a);
407 const int category_b = initial_basis_.GetColumnCategory(col_b);
408 if (category_a != category_b) {
409 return category_a < category_b;
410 } else {
411 return initial_basis_.GetColumnPenalty(col_a) <
412 initial_basis_.GetColumnPenalty(col_b);
413 }
414}
415
416bool InitialBasis::TriangularColumnComparator::operator()(
417 ColIndex col_a, ColIndex col_b) const {
418 if (col_a == col_b) return false;
419 const int category_a = initial_basis_.GetColumnCategory(col_a);
420 const int category_b = initial_basis_.GetColumnCategory(col_b);
421 if (category_a != category_b) {
422 return category_a > category_b;
423 }
424
425 // The nonzero is not in the original Bixby paper, but experiment shows it is
426 // important. It leads to sparser solves, but also sparser direction, which
427 // mean potentially less blocking variables on each pivot...
428 //
429 // TODO(user): Experiments more with this comparator or the
430 // BixbyColumnComparator.
431 if (initial_basis_.compact_matrix_.column(col_a).num_entries() !=
432 initial_basis_.compact_matrix_.column(col_b).num_entries()) {
433 return initial_basis_.compact_matrix_.column(col_a).num_entries() >
434 initial_basis_.compact_matrix_.column(col_b).num_entries();
435 }
436 return initial_basis_.GetColumnPenalty(col_a) >
437 initial_basis_.GetColumnPenalty(col_b);
438}
439
440} // namespace glop
441} // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
#define VLOG(verboselevel)
Definition: base/logging.h:978
ColumnView column(ColIndex col) const
Definition: sparse.h:364
void CompleteTriangularPrimalBasis(ColIndex num_cols, RowToColMapping *basis)
void CompleteTriangularDualBasis(ColIndex num_cols, RowToColMapping *basis)
InitialBasis(const CompactSparseMatrix &compact_matrix, const DenseRow &objective, const DenseRow &lower_bound, const DenseRow &upper_bound, const VariableTypeRow &variable_type)
void CompleteBixbyBasis(ColIndex num_cols, RowToColMapping *basis)
void GetDualMarosBasis(ColIndex num_cols, RowToColMapping *basis)
void GetPrimalMarosBasis(ColIndex num_cols, RowToColMapping *basis)
void ComputeCandidates(ColIndex num_cols, std::vector< ColIndex > *candidates)
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
Fractional InfinityNorm(const DenseColumn &v)
const ColIndex kInvalidCol(-1)
void SetSupportToFalse(const ColumnView &column, DenseBooleanColumn *b)
ColIndex RowToColIndex(RowIndex row)
Definition: lp_types.h:48
bool IsDominated(const ColumnView &column, const DenseColumn &radius)
const RowIndex kInvalidRow(-1)
StrictITIVector< ColIndex, bool > DenseBooleanRow
Definition: lp_types.h:302
StrictITIVector< RowIndex, ColIndex > RowToColMapping
Definition: lp_types.h:342
Fractional RestrictedInfinityNorm(const ColumnView &column, const DenseBooleanColumn &rows_to_consider, RowIndex *row_index)
const double kInfinity
Definition: lp_types.h:83
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
IntVar *const objective_
Definition: search.cc:2951