OR-Tools  8.2
markowitz.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 <limits>
17
18#include "absl/strings/str_format.h"
22
23namespace operations_research {
24namespace glop {
25
27 const CompactSparseMatrixView& basis_matrix, RowPermutation* row_perm,
28 ColumnPermutation* col_perm) {
29 SCOPED_TIME_STAT(&stats_);
30 Clear();
31 const RowIndex num_rows = basis_matrix.num_rows();
32 const ColIndex num_cols = basis_matrix.num_cols();
33 col_perm->assign(num_cols, kInvalidCol);
34 row_perm->assign(num_rows, kInvalidRow);
35
36 // Get the empty matrix corner case out of the way.
37 if (basis_matrix.IsEmpty()) return Status::OK();
38 basis_matrix_ = &basis_matrix;
39
40 // Initialize all the matrices.
41 lower_.Reset(num_rows, num_cols);
42 upper_.Reset(num_rows, num_cols);
43 permuted_lower_.Reset(num_cols);
44 permuted_upper_.Reset(num_cols);
45 permuted_lower_column_needs_solve_.assign(num_cols, false);
46 contains_only_singleton_columns_ = true;
47
48 // Start by moving the singleton columns to the front and by putting their
49 // non-zero coefficient on the diagonal. The general algorithm below would
50 // have the same effect, but this function is a lot faster.
51 int index = 0;
52 ExtractSingletonColumns(basis_matrix, row_perm, col_perm, &index);
53 ExtractResidualSingletonColumns(basis_matrix, row_perm, col_perm, &index);
54 int stats_num_pivots_without_fill_in = index;
55 int stats_degree_two_pivot_columns = 0;
56
57 // Initialize residual_matrix_non_zero_ with the submatrix left after we
58 // removed the singleton and residual singleton columns.
59 residual_matrix_non_zero_.InitializeFromMatrixSubset(
60 basis_matrix, *row_perm, *col_perm, &singleton_column_, &singleton_row_);
61
62 // Perform Gaussian elimination.
63 const int end_index = std::min(num_rows.value(), num_cols.value());
64 const Fractional singularity_threshold =
65 parameters_.markowitz_singularity_threshold();
66 while (index < end_index) {
67 Fractional pivot_coefficient = 0.0;
68 RowIndex pivot_row = kInvalidRow;
69 ColIndex pivot_col = kInvalidCol;
70
71 // TODO(user): If we don't need L and U, we can abort when the residual
72 // matrix becomes dense (i.e. when its density factor is above a certain
73 // threshold). The residual size is 'end_index - index' and the
74 // density can either be computed exactly or estimated from min_markowitz.
75 const int64 min_markowitz = FindPivot(*row_perm, *col_perm, &pivot_row,
76 &pivot_col, &pivot_coefficient);
77
78 // Singular matrix? No pivot will be selected if a column has no entries. If
79 // a column has some entries, then we are sure that a pivot will be selected
80 // but its magnitude can be really close to zero. In both cases, we
81 // report the singularity of the matrix.
82 if (pivot_row == kInvalidRow || pivot_col == kInvalidCol ||
83 std::abs(pivot_coefficient) <= singularity_threshold) {
84 const std::string error_message = absl::StrFormat(
85 "The matrix is singular! pivot = %E", pivot_coefficient);
86 VLOG(1) << "ERROR_LU: " << error_message;
87 return Status(Status::ERROR_LU, error_message);
88 }
89 DCHECK_EQ((*row_perm)[pivot_row], kInvalidRow);
90 DCHECK_EQ((*col_perm)[pivot_col], kInvalidCol);
91
92 // Update residual_matrix_non_zero_.
93 // TODO(user): This step can be skipped, once a fully dense matrix is
94 // obtained. But note that permuted_lower_column_needs_solve_ needs to be
95 // updated.
96 const int pivot_col_degree = residual_matrix_non_zero_.ColDegree(pivot_col);
97 const int pivot_row_degree = residual_matrix_non_zero_.RowDegree(pivot_row);
98 residual_matrix_non_zero_.DeleteRowAndColumn(pivot_row, pivot_col);
99 if (min_markowitz == 0) {
100 ++stats_num_pivots_without_fill_in;
101 if (pivot_col_degree == 1) {
102 RemoveRowFromResidualMatrix(pivot_row, pivot_col);
103 } else {
104 DCHECK_EQ(pivot_row_degree, 1);
105 RemoveColumnFromResidualMatrix(pivot_row, pivot_col);
106 }
107 } else {
108 // TODO(user): Note that in some rare cases, because of numerical
109 // cancellation, the column degree may actually be smaller than
110 // pivot_col_degree. Exploit that better?
112 if (pivot_col_degree == 2) { ++stats_degree_two_pivot_columns; });
113 UpdateResidualMatrix(pivot_row, pivot_col);
114 }
115
116 if (contains_only_singleton_columns_) {
117 DCHECK(permuted_upper_.column(pivot_col).IsEmpty());
118 lower_.AddDiagonalOnlyColumn(1.0);
119 upper_.AddTriangularColumn(basis_matrix.column(pivot_col), pivot_row);
120 } else {
121 lower_.AddAndNormalizeTriangularColumn(permuted_lower_.column(pivot_col),
122 pivot_row, pivot_coefficient);
123 permuted_lower_.ClearAndReleaseColumn(pivot_col);
124
126 permuted_upper_.column(pivot_col), pivot_row, pivot_coefficient);
127 permuted_upper_.ClearAndReleaseColumn(pivot_col);
128 }
129
130 // Update the permutations.
131 (*col_perm)[pivot_col] = ColIndex(index);
132 (*row_perm)[pivot_row] = RowIndex(index);
133 ++index;
134 }
135
136 stats_.pivots_without_fill_in_ratio.Add(
137 1.0 * stats_num_pivots_without_fill_in / end_index);
138 stats_.degree_two_pivot_columns.Add(1.0 * stats_degree_two_pivot_columns /
139 end_index);
140 return Status::OK();
141}
142
144 RowPermutation* row_perm,
145 ColumnPermutation* col_perm,
146 TriangularMatrix* lower, TriangularMatrix* upper) {
147 // The two first swaps allow to use less memory since this way upper_
148 // and lower_ will always stay empty at the end of this function.
149 lower_.Swap(lower);
150 upper_.Swap(upper);
152 ComputeRowAndColumnPermutation(basis_matrix, row_perm, col_perm));
153 SCOPED_TIME_STAT(&stats_);
156 lower_.Swap(lower);
157 upper_.Swap(upper);
158 DCHECK(lower->IsLowerTriangular());
159 DCHECK(upper->IsUpperTriangular());
160 return Status::OK();
161}
162
164 SCOPED_TIME_STAT(&stats_);
165 permuted_lower_.Clear();
166 permuted_upper_.Clear();
167 residual_matrix_non_zero_.Clear();
168 col_by_degree_.Clear();
169 examined_col_.clear();
170 is_col_by_degree_initialized_ = false;
171}
172
173namespace {
174struct MatrixEntry {
175 RowIndex row;
176 ColIndex col;
178 MatrixEntry(RowIndex r, ColIndex c, Fractional coeff)
179 : row(r), col(c), coefficient(coeff) {}
180 bool operator<(const MatrixEntry& o) const {
181 return (row == o.row) ? col < o.col : row < o.row;
182 }
183};
184
185} // namespace
186
187void Markowitz::ExtractSingletonColumns(
188 const CompactSparseMatrixView& basis_matrix, RowPermutation* row_perm,
189 ColumnPermutation* col_perm, int* index) {
190 SCOPED_TIME_STAT(&stats_);
191 std::vector<MatrixEntry> singleton_entries;
192 const ColIndex num_cols = basis_matrix.num_cols();
193 for (ColIndex col(0); col < num_cols; ++col) {
194 const ColumnView& column = basis_matrix.column(col);
195 if (column.num_entries().value() == 1) {
196 singleton_entries.push_back(
197 MatrixEntry(column.GetFirstRow(), col, column.GetFirstCoefficient()));
198 }
199 }
200
201 // Sorting the entries by row indices allows the row_permutation to be closer
202 // to identity which seems like a good idea.
203 std::sort(singleton_entries.begin(), singleton_entries.end());
204 for (const MatrixEntry e : singleton_entries) {
205 if ((*row_perm)[e.row] == kInvalidRow) {
206 (*col_perm)[e.col] = ColIndex(*index);
207 (*row_perm)[e.row] = RowIndex(*index);
208 lower_.AddDiagonalOnlyColumn(1.0);
209 upper_.AddDiagonalOnlyColumn(e.coefficient);
210 ++(*index);
211 }
212 }
213 stats_.basis_singleton_column_ratio.Add(static_cast<double>(*index) /
214 num_cols.value());
215}
216
217bool Markowitz::IsResidualSingletonColumn(const ColumnView& column,
218 const RowPermutation& row_perm,
219 RowIndex* row) {
220 int residual_degree = 0;
221 for (const auto e : column) {
222 if (row_perm[e.row()] != kInvalidRow) continue;
223 ++residual_degree;
224 if (residual_degree > 1) return false;
225 *row = e.row();
226 }
227 return residual_degree == 1;
228}
229
230void Markowitz::ExtractResidualSingletonColumns(
231 const CompactSparseMatrixView& basis_matrix, RowPermutation* row_perm,
232 ColumnPermutation* col_perm, int* index) {
233 SCOPED_TIME_STAT(&stats_);
234 const ColIndex num_cols = basis_matrix.num_cols();
235 RowIndex row = kInvalidRow;
236 for (ColIndex col(0); col < num_cols; ++col) {
237 if ((*col_perm)[col] != kInvalidCol) continue;
238 const ColumnView& column = basis_matrix.column(col);
239 if (!IsResidualSingletonColumn(column, *row_perm, &row)) continue;
240 (*col_perm)[col] = ColIndex(*index);
241 (*row_perm)[row] = RowIndex(*index);
242 lower_.AddDiagonalOnlyColumn(1.0);
243 upper_.AddTriangularColumn(column, row);
244 ++(*index);
245 }
246 stats_.basis_residual_singleton_column_ratio.Add(static_cast<double>(*index) /
247 num_cols.value());
248}
249
250const SparseColumn& Markowitz::ComputeColumn(const RowPermutation& row_perm,
251 ColIndex col) {
252 SCOPED_TIME_STAT(&stats_);
253 // Is this the first time ComputeColumn() sees this column? This is a bit
254 // tricky because just one of the tests is not sufficient in case the matrix
255 // is degenerate.
256 const bool first_time = permuted_lower_.column(col).IsEmpty() &&
257 permuted_upper_.column(col).IsEmpty();
258
259 // If !permuted_lower_column_needs_solve_[col] then the result of the
260 // PermutedLowerSparseSolve() below is already stored in
261 // permuted_lower_.column(col) and we just need to split this column. Note
262 // that this is just an optimization and the code would work if we just
263 // assumed permuted_lower_column_needs_solve_[col] to be always true.
264 SparseColumn* lower_column = permuted_lower_.mutable_column(col);
265 if (permuted_lower_column_needs_solve_[col]) {
266 // Solve a sparse triangular system. If the column 'col' of permuted_lower_
267 // was never computed before by ComputeColumn(), we use the column 'col' of
268 // the matrix to factorize.
269 const ColumnView& input =
270 first_time ? basis_matrix_->column(col) : ColumnView(*lower_column);
271 lower_.PermutedLowerSparseSolve(input, row_perm, lower_column,
272 permuted_upper_.mutable_column(col));
273 permuted_lower_column_needs_solve_[col] = false;
274 return *lower_column;
275 }
276
277 // All the symbolic non-zeros are always present in lower. So if this test is
278 // true, we can conclude that there is no entries from upper that need to be
279 // moved by a cardinality argument.
280 if (lower_column->num_entries() == residual_matrix_non_zero_.ColDegree(col)) {
281 return *lower_column;
282 }
283
284 // In this case, we just need to "split" the lower column. We copy from the
285 // appropriate ColumnView in basis_matrix_.
286 // TODO(user): add PopulateFromColumnView if it is useful elsewhere.
287 if (first_time) {
288 lower_column->Reserve(basis_matrix_->column(col).num_entries());
289 for (const auto e : basis_matrix_->column(col)) {
290 lower_column->SetCoefficient(e.row(), e.coefficient());
291 }
292 }
293 lower_column->MoveTaggedEntriesTo(row_perm,
294 permuted_upper_.mutable_column(col));
295 return *lower_column;
296}
297
298int64 Markowitz::FindPivot(const RowPermutation& row_perm,
299 const ColumnPermutation& col_perm,
300 RowIndex* pivot_row, ColIndex* pivot_col,
301 Fractional* pivot_coefficient) {
302 SCOPED_TIME_STAT(&stats_);
303
304 // Fast track for singleton columns.
305 while (!singleton_column_.empty()) {
306 const ColIndex col = singleton_column_.back();
307 singleton_column_.pop_back();
308 DCHECK_EQ(kInvalidCol, col_perm[col]);
309
310 // This can only happen if the matrix is singular. Continuing will cause
311 // the algorithm to detect the singularity at the end when we stop before
312 // the end.
313 //
314 // TODO(user): We could detect the singularity at this point, but that
315 // may make the code more complex.
316 if (residual_matrix_non_zero_.ColDegree(col) != 1) continue;
317
318 // ComputeColumn() is not used as long as only singleton columns of the
319 // residual matrix are used. See the other condition in
320 // ComputeRowAndColumnPermutation().
321 if (contains_only_singleton_columns_) {
322 *pivot_col = col;
323 for (const SparseColumn::Entry e : basis_matrix_->column(col)) {
324 if (row_perm[e.row()] == kInvalidRow) {
325 *pivot_row = e.row();
326 *pivot_coefficient = e.coefficient();
327 break;
328 }
329 }
330 return 0;
331 }
332 const SparseColumn& column = ComputeColumn(row_perm, col);
333 if (column.IsEmpty()) continue;
334 *pivot_col = col;
335 *pivot_row = column.GetFirstRow();
336 *pivot_coefficient = column.GetFirstCoefficient();
337 return 0;
338 }
339 contains_only_singleton_columns_ = false;
340
341 // Fast track for singleton rows. Note that this is actually more than a fast
342 // track because of the Zlatev heuristic. Such rows may not be processed as
343 // soon as possible otherwise, resulting in more fill-in.
344 while (!singleton_row_.empty()) {
345 const RowIndex row = singleton_row_.back();
346 singleton_row_.pop_back();
347
348 // A singleton row could have been processed when processing a singleton
349 // column. Skip if this is the case.
350 if (row_perm[row] != kInvalidRow) continue;
351
352 // This shows that the matrix is singular, see comment above for the same
353 // case when processing singleton columns.
354 if (residual_matrix_non_zero_.RowDegree(row) != 1) continue;
355 const ColIndex col =
356 residual_matrix_non_zero_.GetFirstNonDeletedColumnFromRow(row);
357 if (col == kInvalidCol) continue;
358 const SparseColumn& column = ComputeColumn(row_perm, col);
359 if (column.IsEmpty()) continue;
360
361 *pivot_col = col;
362 *pivot_row = row;
363 *pivot_coefficient = column.LookUpCoefficient(row);
364 return 0;
365 }
366
367 // col_by_degree_ is not needed before we reach this point. Exploit this with
368 // a lazy initialization.
369 if (!is_col_by_degree_initialized_) {
370 is_col_by_degree_initialized_ = true;
371 const ColIndex num_cols = col_perm.size();
372 col_by_degree_.Reset(row_perm.size().value(), num_cols);
373 for (ColIndex col(0); col < num_cols; ++col) {
374 if (col_perm[col] != kInvalidCol) continue;
375 const int degree = residual_matrix_non_zero_.ColDegree(col);
376 DCHECK_NE(degree, 1);
377 UpdateDegree(col, degree);
378 }
379 }
380
381 // Note(user): we use int64 since this is a product of two ints, moreover
382 // the ints should be relatively small, so that should be fine for a while.
383 int64 min_markowitz_number = std::numeric_limits<int64>::max();
384 examined_col_.clear();
385 const int num_columns_to_examine = parameters_.markowitz_zlatev_parameter();
386 const Fractional threshold = parameters_.lu_factorization_pivot_threshold();
387 while (examined_col_.size() < num_columns_to_examine) {
388 const ColIndex col = col_by_degree_.Pop();
389 if (col == kInvalidCol) break;
390 if (col_perm[col] != kInvalidCol) continue;
391 const int col_degree = residual_matrix_non_zero_.ColDegree(col);
392 examined_col_.push_back(col);
393
394 // Because of the two singleton special cases at the beginning of this
395 // function and because we process columns by increasing degree, we can
396 // derive a lower bound on the best markowitz number we can get by exploring
397 // this column. If we cannot beat this number, we can stop here.
398 //
399 // Note(user): we still process extra column if we can meet the lower bound
400 // to eventually have a better pivot.
401 //
402 // Todo(user): keep the minimum row degree to have a better bound?
403 const int64 markowitz_lower_bound = col_degree - 1;
404 if (min_markowitz_number < markowitz_lower_bound) break;
405
406 // TODO(user): col_degree (which is the same as column.num_entries()) is
407 // actually an upper bound on the number of non-zeros since there may be
408 // numerical cancellations. Exploit this here? Note that it is already used
409 // when we update the non_zero pattern of the residual matrix.
410 const SparseColumn& column = ComputeColumn(row_perm, col);
411 DCHECK_EQ(column.num_entries(), col_degree);
412
413 Fractional max_magnitude = 0.0;
414 for (const SparseColumn::Entry e : column) {
415 max_magnitude = std::max(max_magnitude, std::abs(e.coefficient()));
416 }
417 if (max_magnitude == 0.0) {
418 // All symbolic non-zero entries have been cancelled!
419 // The matrix is singular, but we continue with the other columns.
420 examined_col_.pop_back();
421 continue;
422 }
423
424 const Fractional skip_threshold = threshold * max_magnitude;
425 for (const SparseColumn::Entry e : column) {
426 const Fractional magnitude = std::abs(e.coefficient());
427 if (magnitude < skip_threshold) continue;
428
429 const int row_degree = residual_matrix_non_zero_.RowDegree(e.row());
430 const int64 markowitz_number = (col_degree - 1) * (row_degree - 1);
431 DCHECK_NE(markowitz_number, 0);
432 if (markowitz_number < min_markowitz_number ||
433 ((markowitz_number == min_markowitz_number) &&
434 magnitude > std::abs(*pivot_coefficient))) {
435 min_markowitz_number = markowitz_number;
436 *pivot_col = col;
437 *pivot_row = e.row();
438 *pivot_coefficient = e.coefficient();
439
440 // Note(user): We could abort early here if the markowitz_lower_bound is
441 // reached, but finishing to loop over this column is fast and may lead
442 // to a pivot with a greater magnitude (i.e. a more robust
443 // factorization).
444 }
445 }
446 DCHECK_NE(min_markowitz_number, 0);
447 DCHECK_GE(min_markowitz_number, markowitz_lower_bound);
448 }
449
450 // Push back the columns that we just looked at in the queue since they
451 // are candidates for the next pivot.
452 //
453 // TODO(user): Do that after having updated the matrix? Rationale:
454 // - col_by_degree_ is LIFO, so that may save work in ComputeColumn() by
455 // calling it again on the same columns.
456 // - Maybe the earliest low-degree columns have a better precision? This
457 // actually depends on the number of operations so is not really true.
458 // - Maybe picking the column randomly from the ones with lowest degree would
459 // help having more diversity from one factorization to the next. This is
460 // for the case we do implement this TODO.
461 for (const ColIndex col : examined_col_) {
462 if (col != *pivot_col) {
463 const int degree = residual_matrix_non_zero_.ColDegree(col);
464 col_by_degree_.PushOrAdjust(col, degree);
465 }
466 }
467 return min_markowitz_number;
468}
469
470void Markowitz::UpdateDegree(ColIndex col, int degree) {
471 DCHECK(is_col_by_degree_initialized_);
472
473 // Separating the degree one columns work because we always select such
474 // a column first and pivoting by such columns does not affect the degree of
475 // any other singleton columns (except if the matrix is not inversible).
476 //
477 // Note that using this optimization does change the order in which the
478 // degree one columns are taken compared to pushing them in the queue.
479 if (degree == 1) {
480 // Note that there is no need to remove this column from col_by_degree_
481 // because it will be processed before col_by_degree_.Pop() is called and
482 // then just be ignored.
483 singleton_column_.push_back(col);
484 } else {
485 col_by_degree_.PushOrAdjust(col, degree);
486 }
487}
488
489void Markowitz::RemoveRowFromResidualMatrix(RowIndex pivot_row,
490 ColIndex pivot_col) {
491 SCOPED_TIME_STAT(&stats_);
492 // Note that instead of calling:
493 // residual_matrix_non_zero_.RemoveDeletedColumnsFromRow(pivot_row);
494 // it is a bit faster to test each position with IsColumnDeleted() since we
495 // will not need the pivot row anymore.
496 if (is_col_by_degree_initialized_) {
497 for (const ColIndex col : residual_matrix_non_zero_.RowNonZero(pivot_row)) {
498 if (residual_matrix_non_zero_.IsColumnDeleted(col)) continue;
499 UpdateDegree(col, residual_matrix_non_zero_.DecreaseColDegree(col));
500 }
501 } else {
502 for (const ColIndex col : residual_matrix_non_zero_.RowNonZero(pivot_row)) {
503 if (residual_matrix_non_zero_.IsColumnDeleted(col)) continue;
504 if (residual_matrix_non_zero_.DecreaseColDegree(col) == 1) {
505 singleton_column_.push_back(col);
506 }
507 }
508 }
509}
510
511void Markowitz::RemoveColumnFromResidualMatrix(RowIndex pivot_row,
512 ColIndex pivot_col) {
513 SCOPED_TIME_STAT(&stats_);
514 // The entries of the pivot column are exactly the symbolic non-zeros of the
515 // residual matrix, since we didn't remove the entries with a coefficient of
516 // zero during PermutedLowerSparseSolve().
517 //
518 // Note that it is okay to decrease the degree of a previous pivot row since
519 // it was set to 0 and will never trigger this test. Even if it triggers it,
520 // we just ignore such singleton rows in FindPivot().
521 for (const SparseColumn::Entry e : permuted_lower_.column(pivot_col)) {
522 const RowIndex row = e.row();
523 if (residual_matrix_non_zero_.DecreaseRowDegree(row) == 1) {
524 singleton_row_.push_back(row);
525 }
526 }
527}
528
529void Markowitz::UpdateResidualMatrix(RowIndex pivot_row, ColIndex pivot_col) {
530 SCOPED_TIME_STAT(&stats_);
531 const SparseColumn& pivot_column = permuted_lower_.column(pivot_col);
532 residual_matrix_non_zero_.Update(pivot_row, pivot_col, pivot_column);
533 for (const ColIndex col : residual_matrix_non_zero_.RowNonZero(pivot_row)) {
534 DCHECK_NE(col, pivot_col);
535 UpdateDegree(col, residual_matrix_non_zero_.ColDegree(col));
536 permuted_lower_column_needs_solve_[col] = true;
537 }
538 RemoveColumnFromResidualMatrix(pivot_row, pivot_col);
539}
540
542 row_degree_.clear();
543 col_degree_.clear();
544 row_non_zero_.clear();
545 deleted_columns_.clear();
546 bool_scratchpad_.clear();
547 num_non_deleted_columns_ = 0;
548}
549
550void MatrixNonZeroPattern::Reset(RowIndex num_rows, ColIndex num_cols) {
551 row_degree_.AssignToZero(num_rows);
552 col_degree_.AssignToZero(num_cols);
553 row_non_zero_.clear();
554 row_non_zero_.resize(num_rows.value());
555 deleted_columns_.assign(num_cols, false);
556 bool_scratchpad_.assign(num_cols, false);
557 num_non_deleted_columns_ = num_cols;
558}
559
561 const CompactSparseMatrixView& basis_matrix, const RowPermutation& row_perm,
562 const ColumnPermutation& col_perm, std::vector<ColIndex>* singleton_columns,
563 std::vector<RowIndex>* singleton_rows) {
564 const ColIndex num_cols = basis_matrix.num_cols();
565 const RowIndex num_rows = basis_matrix.num_rows();
566
567 // Reset the matrix and initialize the vectors to the correct sizes.
568 Reset(num_rows, num_cols);
569 singleton_columns->clear();
570 singleton_rows->clear();
571
572 // Compute the number of entries in each row.
573 for (ColIndex col(0); col < num_cols; ++col) {
574 if (col_perm[col] != kInvalidCol) {
575 deleted_columns_[col] = true;
576 --num_non_deleted_columns_;
577 continue;
578 }
579 for (const SparseColumn::Entry e : basis_matrix.column(col)) {
580 ++row_degree_[e.row()];
581 }
582 }
583
584 // Reserve the row_non_zero_ vector sizes.
585 for (RowIndex row(0); row < num_rows; ++row) {
586 if (row_perm[row] == kInvalidRow) {
587 row_non_zero_[row].reserve(row_degree_[row]);
588 if (row_degree_[row] == 1) singleton_rows->push_back(row);
589 } else {
590 // This is needed because in the row degree computation above, we do not
591 // test for row_perm[row] == kInvalidRow because it is a bit faster.
592 row_degree_[row] = 0;
593 }
594 }
595
596 // Initialize row_non_zero_.
597 for (ColIndex col(0); col < num_cols; ++col) {
598 if (col_perm[col] != kInvalidCol) continue;
599 int32 col_degree = 0;
600 for (const SparseColumn::Entry e : basis_matrix.column(col)) {
601 const RowIndex row = e.row();
602 if (row_perm[row] == kInvalidRow) {
603 ++col_degree;
604 row_non_zero_[row].push_back(col);
605 }
606 }
607 col_degree_[col] = col_degree;
608 if (col_degree == 1) singleton_columns->push_back(col);
609 }
610}
611
612void MatrixNonZeroPattern::AddEntry(RowIndex row, ColIndex col) {
613 ++row_degree_[row];
614 ++col_degree_[col];
615 row_non_zero_[row].push_back(col);
616}
617
619 return --col_degree_[col];
620}
621
623 return --row_degree_[row];
624}
625
627 ColIndex pivot_col) {
628 DCHECK(!deleted_columns_[pivot_col]);
629 deleted_columns_[pivot_col] = true;
630 --num_non_deleted_columns_;
631
632 // We do that to optimize RemoveColumnFromResidualMatrix().
633 row_degree_[pivot_row] = 0;
634}
635
637 return deleted_columns_[col];
638}
639
641 auto& ref = row_non_zero_[row];
642 int new_index = 0;
643 const int end = ref.size();
644 for (int i = 0; i < end; ++i) {
645 const ColIndex col = ref[i];
646 if (!deleted_columns_[col]) {
647 ref[new_index] = col;
648 ++new_index;
649 }
650 }
651 ref.resize(new_index);
652}
653
655 RowIndex row) const {
656 for (const ColIndex col : RowNonZero(row)) {
657 if (!IsColumnDeleted(col)) return col;
658 }
659 return kInvalidCol;
660}
661
662void MatrixNonZeroPattern::Update(RowIndex pivot_row, ColIndex pivot_col,
663 const SparseColumn& column) {
664 // Since DeleteRowAndColumn() must be called just before this function,
665 // the pivot column has been marked as deleted but degrees have not been
666 // updated yet. Hence the +1.
667 DCHECK(deleted_columns_[pivot_col]);
668 const int max_row_degree = num_non_deleted_columns_.value() + 1;
669
671 for (const ColIndex col : row_non_zero_[pivot_row]) {
673 bool_scratchpad_[col] = false;
674 }
675
676 // We only need to merge the row for the position with a coefficient different
677 // from 0.0. Note that the column must contain all the symbolic non-zeros for
678 // the row degree to be updated correctly. Note also that decreasing the row
679 // degrees due to the deletion of pivot_col will happen outside this function.
680 for (const SparseColumn::Entry e : column) {
681 const RowIndex row = e.row();
682 if (row == pivot_row) continue;
683
684 // If the row is fully dense, there is nothing to do (the merge below will
685 // not change anything). This is a small price to pay for a huge gain when
686 // the matrix becomes dense.
687 if (e.coefficient() == 0.0 || row_degree_[row] == max_row_degree) continue;
688 DCHECK_LT(row_degree_[row], max_row_degree);
689
690 // We only clean row_non_zero_[row] if there are more than 4 entries to
691 // delete. Note(user): the 4 is somewhat arbitrary, but gives good results
692 // on the Netlib (23/04/2013). Note that calling
693 // RemoveDeletedColumnsFromRow() is not mandatory and does not change the LU
694 // decomposition, so we could call it all the time or never and the
695 // algorithm would still work.
696 const int kDeletionThreshold = 4;
697 if (row_non_zero_[row].size() > row_degree_[row] + kDeletionThreshold) {
699 }
700 // TODO(user): Special case if row_non_zero_[pivot_row].size() == 1?
701 if (/* DISABLES CODE */ (true)) {
702 MergeInto(pivot_row, row);
703 } else {
704 // This is currently not used, but kept as an alternative algorithm to
705 // investigate. The performance is really similar, but the final L.U is
706 // different. Note that when this is used, there is no need to modify
707 // bool_scratchpad_ at the beginning of this function.
708 //
709 // TODO(user): Add unit tests before using this.
710 MergeIntoSorted(pivot_row, row);
711 }
712 }
713}
714
715void MatrixNonZeroPattern::MergeInto(RowIndex pivot_row, RowIndex row) {
716 // Note that bool_scratchpad_ must be already false on the positions in
717 // row_non_zero_[pivot_row].
718 for (const ColIndex col : row_non_zero_[row]) {
719 bool_scratchpad_[col] = true;
720 }
721
722 auto& non_zero = row_non_zero_[row];
723 const int old_size = non_zero.size();
724 for (const ColIndex col : row_non_zero_[pivot_row]) {
725 if (bool_scratchpad_[col]) {
726 bool_scratchpad_[col] = false;
727 } else {
728 non_zero.push_back(col);
729 ++col_degree_[col];
730 }
731 }
732 row_degree_[row] += non_zero.size() - old_size;
733}
734
735namespace {
736
737// Given two sorted vectors (the second one is the initial value of out), merges
738// them and outputs the sorted result in out. The merge is stable and an element
739// of input_a will appear before the identical elements of the second input.
740template <typename V, typename W>
741void MergeSortedVectors(const V& input_a, W* out) {
742 if (input_a.empty()) return;
743 const auto& input_b = *out;
744 int index_a = input_a.size() - 1;
745 int index_b = input_b.size() - 1;
746 int index_out = input_a.size() + input_b.size();
747 out->resize(index_out);
748 while (index_a >= 0) {
749 if (index_b < 0) {
750 while (index_a >= 0) {
751 --index_out;
752 (*out)[index_out] = input_a[index_a];
753 --index_a;
754 }
755 return;
756 }
757 --index_out;
758 if (input_a[index_a] > input_b[index_b]) {
759 (*out)[index_out] = input_a[index_a];
760 --index_a;
761 } else {
762 (*out)[index_out] = input_b[index_b];
763 --index_b;
764 }
765 }
766}
767
768} // namespace
769
770// The algorithm first computes into col_scratchpad_ the entries in pivot_row
771// that are not in the row (i.e. the fill-in). It then updates the non-zero
772// pattern using this temporary vector.
773void MatrixNonZeroPattern::MergeIntoSorted(RowIndex pivot_row, RowIndex row) {
774 // We want to add the entries of the input not already in the output.
775 const auto& input = row_non_zero_[pivot_row];
776 const auto& output = row_non_zero_[row];
777
778 // These two resizes are because of the set_difference() output iterator api.
779 col_scratchpad_.resize(input.size());
780 col_scratchpad_.resize(std::set_difference(input.begin(), input.end(),
781 output.begin(), output.end(),
782 col_scratchpad_.begin()) -
783 col_scratchpad_.begin());
784
785 // Add the fill-in to the pattern.
786 for (const ColIndex col : col_scratchpad_) {
787 ++col_degree_[col];
788 }
789 row_degree_[row] += col_scratchpad_.size();
790 MergeSortedVectors(col_scratchpad_, &row_non_zero_[row]);
791}
792
794 col_degree_.clear();
795 col_index_.clear();
796 col_by_degree_.clear();
797}
798
799void ColumnPriorityQueue::Reset(int max_degree, ColIndex num_cols) {
800 Clear();
801 col_degree_.assign(num_cols, 0);
802 col_index_.assign(num_cols, -1);
803 col_by_degree_.resize(max_degree + 1);
804 min_degree_ = num_cols.value();
805}
806
808 DCHECK_GE(degree, 0);
809 DCHECK_LT(degree, col_by_degree_.size());
810 const int32 old_degree = col_degree_[col];
811 if (degree != old_degree) {
812 const int32 old_index = col_index_[col];
813 if (old_index != -1) {
814 col_by_degree_[old_degree][old_index] = col_by_degree_[old_degree].back();
815 col_index_[col_by_degree_[old_degree].back()] = old_index;
816 col_by_degree_[old_degree].pop_back();
817 }
818 if (degree > 0) {
819 col_index_[col] = col_by_degree_[degree].size();
820 col_degree_[col] = degree;
821 col_by_degree_[degree].push_back(col);
822 min_degree_ = std::min(min_degree_, degree);
823 } else {
824 col_index_[col] = -1;
825 col_degree_[col] = 0;
826 }
827 }
828}
829
831 while (col_by_degree_[min_degree_].empty()) {
832 min_degree_++;
833 if (min_degree_ == col_by_degree_.size()) return kInvalidCol;
834 }
835 const ColIndex col = col_by_degree_[min_degree_].back();
836 col_by_degree_[min_degree_].pop_back();
837 col_index_[col] = -1;
838 col_degree_[col] = 0;
839 return col;
840}
841
843 mapping_.assign(num_cols.value(), -1);
844 free_columns_.clear();
845 columns_.clear();
846}
847
849 ColIndex col) const {
850 if (mapping_[col] == -1) return empty_column_;
851 return columns_[mapping_[col]];
852}
853
855 ColIndex col) {
856 if (mapping_[col] != -1) return &columns_[mapping_[col]];
857 int new_col_index;
858 if (free_columns_.empty()) {
859 new_col_index = columns_.size();
860 columns_.push_back(SparseColumn());
861 } else {
862 new_col_index = free_columns_.back();
863 free_columns_.pop_back();
864 }
865 mapping_[col] = new_col_index;
866 return &columns_[new_col_index];
867}
868
870 DCHECK_NE(mapping_[col], -1);
871 free_columns_.push_back(mapping_[col]);
872 columns_[mapping_[col]].Clear();
873 mapping_[col] = -1;
874}
875
877 mapping_.clear();
878 free_columns_.clear();
879 columns_.clear();
880}
881
882} // namespace glop
883} // 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_GE(val1, val2)
Definition: base/logging.h:889
#define DCHECK_LT(val1, val2)
Definition: base/logging.h:888
#define DCHECK(condition)
Definition: base/logging.h:884
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
#define VLOG(verboselevel)
Definition: base/logging.h:978
void assign(size_type n, const value_type &val)
void resize(size_type new_size)
void reserve(size_type n)
size_type size() const
void push_back(const value_type &x)
void Reset(int32 max_degree, ColIndex num_cols)
Definition: markowitz.cc:799
void PushOrAdjust(ColIndex col, int32 degree)
Definition: markowitz.cc:807
const ColumnView column(ColIndex col) const
Definition: sparse.h:481
ABSL_MUST_USE_RESULT Status ComputeLU(const CompactSparseMatrixView &basis_matrix, RowPermutation *row_perm, ColumnPermutation *col_perm, TriangularMatrix *lower, TriangularMatrix *upper)
Definition: markowitz.cc:143
ABSL_MUST_USE_RESULT Status ComputeRowAndColumnPermutation(const CompactSparseMatrixView &basis_matrix, RowPermutation *row_perm, ColumnPermutation *col_perm)
Definition: markowitz.cc:26
const absl::InlinedVector< ColIndex, 6 > & RowNonZero(RowIndex row) const
Definition: markowitz.h:167
void DeleteRowAndColumn(RowIndex pivot_row, ColIndex pivot_col)
Definition: markowitz.cc:626
void AddEntry(RowIndex row, ColIndex col)
Definition: markowitz.cc:612
void Reset(RowIndex num_rows, ColIndex num_cols)
Definition: markowitz.cc:550
void Update(RowIndex pivot_row, ColIndex pivot_col, const SparseColumn &column)
Definition: markowitz.cc:662
ColIndex GetFirstNonDeletedColumnFromRow(RowIndex row) const
Definition: markowitz.cc:654
void InitializeFromMatrixSubset(const CompactSparseMatrixView &basis_matrix, const RowPermutation &row_perm, const ColumnPermutation &col_perm, std::vector< ColIndex > *singleton_columns, std::vector< RowIndex > *singleton_rows)
Definition: markowitz.cc:560
void assign(IndexType size, IndexType value)
const SparseColumn & column(ColIndex col) const
Definition: markowitz.cc:848
static const Status OK()
Definition: status.h:54
void assign(IntType size, const T &v)
Definition: lp_types.h:274
void AddTriangularColumnWithGivenDiagonalEntry(const SparseColumn &column, RowIndex diagonal_row, Fractional diagonal_value)
Definition: sparse.cc:672
void Swap(TriangularMatrix *other)
Definition: sparse.cc:593
void AddAndNormalizeTriangularColumn(const SparseColumn &column, RowIndex diagonal_row, Fractional diagonal_coefficient)
Definition: sparse.cc:655
void AddTriangularColumn(const ColumnView &column, RowIndex diagonal_row)
Definition: sparse.cc:640
void PermutedLowerSparseSolve(const ColumnView &rhs, const RowPermutation &row_perm, SparseColumn *lower, SparseColumn *upper)
Definition: sparse.cc:1038
void ApplyRowPermutationToNonDiagonalEntries(const RowPermutation &row_perm)
Definition: sparse.cc:712
void AddDiagonalOnlyColumn(Fractional diagonal_value)
Definition: sparse.cc:636
void Reset(RowIndex num_rows, ColIndex col_capacity)
Definition: sparse.cc:524
int int32
int64_t int64
Fractional coefficient
Definition: markowitz.cc:177
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
const ColIndex kInvalidCol(-1)
Permutation< ColIndex > ColumnPermutation
const RowIndex kInvalidRow(-1)
Permutation< RowIndex > RowPermutation
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int index
Definition: pack.cc:508
static int input(yyscan_t yyscanner)
#define IF_STATS_ENABLED(instructions)
Definition: stats.h:437
#define SCOPED_TIME_STAT(stats)
Definition: stats.h:438
#define GLOP_RETURN_IF_ERROR(function_call)
Definition: status.h:70