OR-Tools  8.1
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 
14 #include "ortools/glop/markowitz.h"
15 
16 #include <limits>
17 
18 #include "absl/strings/str_format.h"
21 #include "ortools/lp_data/sparse.h"
22 
23 namespace operations_research {
24 namespace 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_);
154  lower_.ApplyRowPermutationToNonDiagonalEntries(*row_perm);
155  upper_.ApplyRowPermutationToNonDiagonalEntries(*row_perm);
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 
173 namespace {
174 struct 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 
187 void 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 
217 bool 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 
230 void 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 
250 const 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 
298 int64 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 
470 void 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 
489 void 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 
511 void 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 
529 void 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 
550 void 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 
612 void 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 
662 void 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 
670  RemoveDeletedColumnsFromRow(pivot_row);
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 
715 void 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 
735 namespace {
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.
740 template <typename V, typename W>
741 void 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.
773 void 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 
799 void 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
operations_research::glop::MatrixNonZeroPattern::RowDegree
int32 RowDegree(RowIndex row) const
Definition: markowitz.h:161
operations_research::glop::MatrixNonZeroPattern::Reset
void Reset(RowIndex num_rows, ColIndex num_cols)
Definition: markowitz.cc:550
min
int64 min
Definition: alldiff_cst.cc:138
DCHECK_LT
#define DCHECK_LT(val1, val2)
Definition: base/logging.h:888
VLOG
#define VLOG(verboselevel)
Definition: base/logging.h:978
operations_research::glop::TriangularMatrix::Swap
void Swap(TriangularMatrix *other)
Definition: sparse.cc:593
absl::StrongVector::push_back
void push_back(const value_type &x)
Definition: strong_vector.h:158
operations_research::glop::MatrixNonZeroPattern::Clear
void Clear()
Definition: markowitz.cc:541
max
int64 max
Definition: alldiff_cst.cc:139
operations_research::glop::TriangularMatrix::AddTriangularColumnWithGivenDiagonalEntry
void AddTriangularColumnWithGivenDiagonalEntry(const SparseColumn &column, RowIndex diagonal_row, Fractional diagonal_value)
Definition: sparse.cc:672
IF_STATS_ENABLED
#define IF_STATS_ENABLED(instructions)
Definition: stats.h:435
operations_research::glop::kInvalidCol
const ColIndex kInvalidCol(-1)
absl::StrongVector::size
size_type size() const
Definition: strong_vector.h:147
operations_research::glop::SparseMatrixWithReusableColumnMemory::mutable_column
SparseColumn * mutable_column(ColIndex col)
Definition: markowitz.cc:854
operations_research::glop::MatrixNonZeroPattern::RemoveDeletedColumnsFromRow
void RemoveDeletedColumnsFromRow(RowIndex row)
Definition: markowitz.cc:640
operations_research::glop::TriangularMatrix::AddAndNormalizeTriangularColumn
void AddAndNormalizeTriangularColumn(const SparseColumn &column, RowIndex diagonal_row, Fractional diagonal_coefficient)
Definition: sparse.cc:655
operations_research::glop::CompactSparseMatrixView
Definition: sparse.h:471
operations_research::glop::MatrixNonZeroPattern::Update
void Update(RowIndex pivot_row, ColIndex pivot_col, const SparseColumn &column)
Definition: markowitz.cc:662
operations_research::glop::Status
Definition: status.h:24
operations_research::glop::ColumnPriorityQueue::Pop
ColIndex Pop()
Definition: markowitz.cc:830
lp_utils.h
operations_research::glop::MatrixNonZeroPattern::DeleteRowAndColumn
void DeleteRowAndColumn(RowIndex pivot_row, ColIndex pivot_col)
Definition: markowitz.cc:626
operations_research::glop::StrictITIVector::size
IntType size() const
Definition: lp_types.h:276
operations_research
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
Definition: dense_doubly_linked_list.h:21
operations_research::glop::TriangularMatrix::AddTriangularColumn
void AddTriangularColumn(const ColumnView &column, RowIndex diagonal_row)
Definition: sparse.cc:640
operations_research::glop::Status::OK
static const Status OK()
Definition: status.h:54
operations_research::glop::SparseColumn
Definition: sparse_column.h:44
int64
int64_t int64
Definition: integral_types.h:34
index
int index
Definition: pack.cc:508
int32
int int32
Definition: integral_types.h:33
operations_research::glop::SparseMatrixWithReusableColumnMemory::Reset
void Reset(ColIndex num_cols)
Definition: markowitz.cc:842
operations_research::glop::TriangularMatrix::AddDiagonalOnlyColumn
void AddDiagonalOnlyColumn(Fractional diagonal_value)
Definition: sparse.cc:636
absl::StrongVector::assign
void assign(size_type n, const value_type &val)
Definition: strong_vector.h:131
operations_research::glop::Fractional
double Fractional
Definition: lp_types.h:77
operations_research::glop::MatrixNonZeroPattern::IsColumnDeleted
bool IsColumnDeleted(ColIndex col) const
Definition: markowitz.cc:636
DCHECK_NE
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
operations_research::glop::CompactSparseMatrixView::column
const ColumnView column(ColIndex col) const
Definition: sparse.h:481
operations_research::glop::ColumnPriorityQueue::PushOrAdjust
void PushOrAdjust(ColIndex col, int32 degree)
Definition: markowitz.cc:807
SCOPED_TIME_STAT
#define SCOPED_TIME_STAT(stats)
Definition: stats.h:436
operations_research::glop::TriangularMatrix::Reset
void Reset(RowIndex num_rows, ColIndex col_capacity)
Definition: sparse.cc:524
operations_research::glop::Status::ERROR_LU
@ ERROR_LU
Definition: status.h:32
operations_research::glop::TriangularMatrix::PermutedLowerSparseSolve
void PermutedLowerSparseSolve(const ColumnView &rhs, const RowPermutation &row_perm, SparseColumn *lower, SparseColumn *upper)
Definition: sparse.cc:1038
operations_research::glop::StrictITIVector::assign
void assign(IntType size, const T &v)
Definition: lp_types.h:274
absl::StrongVector::reserve
void reserve(size_type n)
Definition: strong_vector.h:157
operations_research::glop::MatrixNonZeroPattern::GetFirstNonDeletedColumnFromRow
ColIndex GetFirstNonDeletedColumnFromRow(RowIndex row) const
Definition: markowitz.cc:654
operations_research::glop::SparseMatrixWithReusableColumnMemory::Clear
void Clear()
Definition: markowitz.cc:876
absl::StrongVector::pop_back
void pop_back()
Definition: strong_vector.h:168
operations_research::glop::Markowitz::ComputeRowAndColumnPermutation
ABSL_MUST_USE_RESULT Status ComputeRowAndColumnPermutation(const CompactSparseMatrixView &basis_matrix, RowPermutation *row_perm, ColumnPermutation *col_perm)
Definition: markowitz.cc:26
operations_research::glop::MatrixNonZeroPattern::InitializeFromMatrixSubset
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
operations_research::glop::CompactSparseMatrixView::num_rows
RowIndex num_rows() const
Definition: sparse.h:479
operations_research::glop::MatrixNonZeroPattern::ColDegree
int32 ColDegree(ColIndex col) const
Definition: markowitz.h:154
operations_research::glop::MatrixNonZeroPattern::AddEntry
void AddEntry(RowIndex row, ColIndex col)
Definition: markowitz.cc:612
operations_research::glop::ColumnPriorityQueue::Reset
void Reset(int32 max_degree, ColIndex num_cols)
Definition: markowitz.cc:799
operations_research::glop::Permutation::assign
void assign(IndexType size, IndexType value)
Definition: lp_data/permutation.h:58
operations_research::glop::StrictITIVector::AssignToZero
void AssignToZero(IntType size)
Definition: lp_types.h:290
DCHECK
#define DCHECK(condition)
Definition: base/logging.h:884
operations_research::glop::MatrixNonZeroPattern::DecreaseRowDegree
int32 DecreaseRowDegree(RowIndex row)
Definition: markowitz.cc:622
operations_research::glop::TriangularMatrix
Definition: sparse.h:502
operations_research::glop::RowPermutation
Permutation< RowIndex > RowPermutation
Definition: lp_data/permutation.h:93
markowitz.h
operations_research::glop::Permutation< RowIndex >
operations_research::glop::TriangularMatrix::IsUpperTriangular
bool IsUpperTriangular() const
Definition: sparse.cc:702
operations_research::glop::Markowitz::Clear
void Clear()
Definition: markowitz.cc:163
operations_research::glop::ColumnPermutation
Permutation< ColIndex > ColumnPermutation
Definition: lp_data/permutation.h:94
operations_research::glop::TriangularMatrix::IsLowerTriangular
bool IsLowerTriangular() const
Definition: sparse.cc:692
operations_research::glop::SparseVector< RowIndex, SparseColumnIterator >::Entry
typename Iterator::Entry Entry
Definition: sparse_vector.h:91
DCHECK_GE
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
operations_research::glop::ColumnView::num_entries
EntryIndex num_entries() const
Definition: sparse_column.h:82
absl::StrongVector::clear
void clear()
Definition: strong_vector.h:170
operations_research::glop::CompactSparseMatrixView::IsEmpty
bool IsEmpty() const
Definition: sparse.h:478
operations_research::glop::SparseMatrixWithReusableColumnMemory::column
const SparseColumn & column(ColIndex col) const
Definition: markowitz.cc:848
col
ColIndex col
Definition: markowitz.cc:176
operations_research::glop::SparseVector::IsEmpty
bool IsEmpty() const
Definition: sparse_vector.h:529
operations_research::glop::SparseMatrixWithReusableColumnMemory::ClearAndReleaseColumn
void ClearAndReleaseColumn(ColIndex col)
Definition: markowitz.cc:869
input
static int input(yyscan_t yyscanner)
row
RowIndex row
Definition: markowitz.cc:175
DCHECK_EQ
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
operations_research::glop::ColumnPriorityQueue::Clear
void Clear()
Definition: markowitz.cc:793
absl::StrongVector::resize
void resize(size_type new_size)
Definition: strong_vector.h:150
operations_research::glop::TriangularMatrix::ApplyRowPermutationToNonDiagonalEntries
void ApplyRowPermutationToNonDiagonalEntries(const RowPermutation &row_perm)
Definition: sparse.cc:712
operations_research::glop::MatrixNonZeroPattern::DecreaseColDegree
int32 DecreaseColDegree(ColIndex col)
Definition: markowitz.cc:618
lp_types.h
operations_research::glop::kInvalidRow
const RowIndex kInvalidRow(-1)
coefficient
Fractional coefficient
Definition: markowitz.cc:177
operations_research::glop::Markowitz::ComputeLU
ABSL_MUST_USE_RESULT Status ComputeLU(const CompactSparseMatrixView &basis_matrix, RowPermutation *row_perm, ColumnPermutation *col_perm, TriangularMatrix *lower, TriangularMatrix *upper)
Definition: markowitz.cc:143
operations_research::glop::CompactSparseMatrixView::num_cols
ColIndex num_cols() const
Definition: sparse.h:480
sparse.h
GLOP_RETURN_IF_ERROR
#define GLOP_RETURN_IF_ERROR(function_call)
Definition: status.h:70
operations_research::glop::MatrixNonZeroPattern::RowNonZero
const absl::InlinedVector< ColIndex, 6 > & RowNonZero(RowIndex row) const
Definition: markowitz.h:166