OR-Tools  8.2
markowitz.h
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 // LU decomposition algorithm of a sparse matrix B with Markowitz pivot
15 // selection strategy. The algorithm constructs a lower matrix L, upper matrix
16 // U, row permutation P and a column permutation Q such that L.U = P.B.Q^{-1}.
17 //
18 // The current algorithm is a mix of ideas that can be found in the literature
19 // and of some optimizations tailored for its use in a revised simplex algorithm
20 // (like a fast processing of the singleton columns present in B). It constructs
21 // L and U column by column from left to right.
22 //
23 // A key concept is the one of the residual matrix which is the bottom right
24 // square submatrix that still needs to be factorized during the classical
25 // Gaussian elimination. The algorithm maintains the non-zero pattern of its
26 // rows and its row/column degrees.
27 //
28 // At each step, a number of columns equal to 'markowitz_zlatev_parameter' are
29 // chosen as candidates from the residual matrix. They are the ones with minimal
30 // residual column degree. They can be found easily because the columns of the
31 // residual matrix are kept in a priority queue.
32 //
33 // We compute the numerical value of these residual columns like in a
34 // left-looking algorithm by solving a sparse lower-triangular system with the
35 // current L constructed so far. Note that this step is highly optimized for
36 // sparsity and we reuse the computations done in the previous steps (if the
37 // candidate column was already considered before). As a by-product, we also
38 // get the corresponding column of U.
39 //
40 // Among the entries of these columns, a pivot is chosen such that the product:
41 // (num_column_entries - 1) * (num_row_entries - 1)
42 // is minimized. Only the pivots with a magnitude greater than
43 // 'lu_factorization_pivot_threshold' times the maximum magnitude of the
44 // corresponding residual column are considered for stability reasons.
45 //
46 // Once the pivot is chosen, the residual column divided by the pivot becomes a
47 // column of L, and the non-zero pattern of the new residual submatrix is
48 // updated by subtracting the outer product of this pivot column times the pivot
49 // row. The product minimized above is thus an upper bound of the number of
50 // fill-in created during a step.
51 //
52 // References:
53 //
54 // J. R. Gilbert and T. Peierls, "Sparse partial pivoting in time proportional
55 // to arithmetic operations," SIAM J. Sci. Statist. Comput., 9 (1988): 862-874.
56 //
57 // I.S. Duff, A.M. Erisman and J.K. Reid, "Direct Methods for Sparse Matrices",
58 // Clarendon, Oxford, UK, 1987, ISBN 0-19-853421-3,
59 // http://www.amazon.com/dp/0198534213
60 //
61 // T.A. Davis, "Direct methods for Sparse Linear Systems", SIAM, Philadelphia,
62 // 2006, ISBN-13: 978-0-898716-13, http://www.amazon.com/dp/0898716136
63 //
64 // TODO(user): Determine whether any of these would bring any benefit:
65 // - S.C. Eisenstat and J.W.H. Liu, "The theory of elimination trees for
66 // sparse unsymmetric matrices," SIAM J. Matrix Anal. Appl., 26:686-705,
67 // January 2005
68 // - S.C. Eisenstat and J.W.H. Liu. "Algorithmic aspects of elimination trees
69 // for sparse unsymmetric matrices," SIAM J. Matrix Anal. Appl.,
70 // 29:1363-1381, January 2008.
71 // - http://perso.ens-lyon.fr/~bucar/papers/kauc.pdf
72 
73 #ifndef OR_TOOLS_GLOP_MARKOWITZ_H_
74 #define OR_TOOLS_GLOP_MARKOWITZ_H_
75 
76 #include <queue>
77 
78 #include "absl/container/inlined_vector.h"
79 #include "ortools/base/logging.h"
81 #include "ortools/glop/parameters.pb.h"
82 #include "ortools/glop/status.h"
85 #include "ortools/lp_data/sparse.h"
87 #include "ortools/util/stats.h"
88 
89 namespace operations_research {
90 namespace glop {
91 
92 // Holds the non-zero positions (by row) and column/row degree of the residual
93 // matrix during the Gaussian elimination.
94 //
95 // During each step of Gaussian elimination, a row and a column will be
96 // "removed" from the residual matrix. Note however that the row and column
97 // indices of the non-removed part do not change, so the residual matrix at a
98 // given step will only correspond to a subset of the initial indices.
100  public:
102 
103  // Releases the memory used by this class.
104  void Clear();
105 
106  // Resets the pattern to the one of an empty square matrix of the given size.
107  void Reset(RowIndex num_rows, ColIndex num_cols);
108 
109  // Resets the pattern to the one of the given matrix but only for the
110  // rows/columns whose given permutation is kInvalidRow or kInvalidCol.
111  // This also fills the singleton columns/rows with the corresponding entries.
112  void InitializeFromMatrixSubset(const CompactSparseMatrixView& basis_matrix,
113  const RowPermutation& row_perm,
114  const ColumnPermutation& col_perm,
115  std::vector<ColIndex>* singleton_columns,
116  std::vector<RowIndex>* singleton_rows);
117 
118  // Adds a non-zero entry to the matrix. There should be no duplicates.
119  void AddEntry(RowIndex row, ColIndex col);
120 
121  // Marks the given pivot row and column as deleted.
122  // This is called at each step of the Gaussian elimination on the pivot.
123  void DeleteRowAndColumn(RowIndex pivot_row, ColIndex pivot_col);
124 
125  // Decreases the degree of a row/column. This is the basic operation used to
126  // keep the correct degree after a call to DeleteRowAndColumn(). This is
127  // because row_non_zero_[row] is only lazily cleaned.
128  int32 DecreaseRowDegree(RowIndex row);
129  int32 DecreaseColDegree(ColIndex col);
130 
131  // Returns true if the column has been deleted by DeleteRowAndColumn().
132  bool IsColumnDeleted(ColIndex col) const;
133 
134  // Removes from the corresponding row_non_zero_[row] the columns that have
135  // been previously deleted by DeleteRowAndColumn().
136  void RemoveDeletedColumnsFromRow(RowIndex row);
137 
138  // Returns the first non-deleted column index from this row or kInvalidCol if
139  // none can be found.
140  ColIndex GetFirstNonDeletedColumnFromRow(RowIndex row) const;
141 
142  // Performs a generic Gaussian update of the residual matrix:
143  // - DeleteRowAndColumn() must already have been called.
144  // - The non-zero pattern is augmented (set union) by the one of the
145  // outer product of the pivot column and row.
146  //
147  // Important: as a small optimization, this function does not call
148  // DecreaseRowDegree() on the row in the pivot column. This has to be done by
149  // the client.
150  void Update(RowIndex pivot_row, ColIndex pivot_col,
151  const SparseColumn& column);
152 
153  // Returns the degree (i.e. the number of non-zeros) of the given column.
154  // This is only valid for the column indices still in the residual matrix.
155  int32 ColDegree(ColIndex col) const {
156  DCHECK(!deleted_columns_[col]);
157  return col_degree_[col];
158  }
159 
160  // Returns the degree (i.e. the number of non-zeros) of the given row.
161  // This is only valid for the row indices still in the residual matrix.
162  int32 RowDegree(RowIndex row) const { return row_degree_[row]; }
163 
164  // Returns the set of non-zeros of the given row (unsorted).
165  // Call RemoveDeletedColumnsFromRow(row) to clean the row first.
166  // This is only valid for the row indices still in the residual matrix.
167  const absl::InlinedVector<ColIndex, 6>& RowNonZero(RowIndex row) const {
168  return row_non_zero_[row];
169  }
170 
171  private:
172  // Augments the non-zero pattern of the given row by taking its union with the
173  // non-zero pattern of the given pivot_row.
174  void MergeInto(RowIndex pivot_row, RowIndex row);
175 
176  // Different version of MergeInto() that works only if the non-zeros position
177  // of each row are sorted in increasing order. The output will also be sorted.
178  //
179  // TODO(user): This is currently not used but about the same speed as the
180  // non-sorted version. Investigate more.
181  void MergeIntoSorted(RowIndex pivot_row, RowIndex row);
182 
183  // Using InlinedVector helps because we usually have many rows with just a few
184  // non-zeros. Note that on a 64 bits computer we get exactly 6 inlined int32
185  // elements without extra space, and the size of the inlined vector is 4 times
186  // 64 bits.
187  //
188  // TODO(user): We could be even more efficient since a size of int32 is enough
189  // for us and we could store in common the inlined/not-inlined size.
193  DenseBooleanRow deleted_columns_;
194  DenseBooleanRow bool_scratchpad_;
195  std::vector<ColIndex> col_scratchpad_;
196  ColIndex num_non_deleted_columns_;
197 
198  DISALLOW_COPY_AND_ASSIGN(MatrixNonZeroPattern);
199 };
200 
201 // Adjustable priority queue of columns. Pop() returns a column with the
202 // smallest degree first (degree = number of entries in the column).
203 // Empty columns (i.e. with degree 0) are not stored in the queue.
205  public:
207 
208  // Releases the memory used by this class.
209  void Clear();
210 
211  // Clears the queue and prepares it to store up to num_cols column indices
212  // with a degree from 1 to max_degree included.
213  void Reset(int32 max_degree, ColIndex num_cols);
214 
215  // Changes the degree of a column and make sure it is in the queue. The degree
216  // must be non-negative (>= 0) and at most equal to the value of num_cols used
217  // in Reset(). A degree of zero will remove the column from the queue.
218  void PushOrAdjust(ColIndex col, int32 degree);
219 
220  // Removes the column index with higher priority from the queue and returns
221  // it. Returns kInvalidCol if the queue is empty.
222  ColIndex Pop();
223 
224  private:
227  std::vector<std::vector<ColIndex>> col_by_degree_;
228  int32 min_degree_;
229  DISALLOW_COPY_AND_ASSIGN(ColumnPriorityQueue);
230 };
231 
232 // Contains a set of columns indexed by ColIndex. This is like a SparseMatrix
233 // but this class is optimized for the case where only a small subset of columns
234 // is needed at the same time (like it is the case in our LU algorithm). It
235 // reuses the memory of the columns that are no longer needed.
237  public:
239 
240  // Resets the repository to num_cols empty columns.
241  void Reset(ColIndex num_cols);
242 
243  // Returns the column with given index.
244  const SparseColumn& column(ColIndex col) const;
245 
246  // Gets the mutable column with given column index. The returned vector
247  // address is only valid until the next call to mutable_column().
248  SparseColumn* mutable_column(ColIndex col);
249 
250  // Clears the column with given index and releases its memory to the common
251  // memory pool that is used to create new mutable_column() on demand.
252  void ClearAndReleaseColumn(ColIndex col);
253 
254  // Reverts this class to its initial state. This releases the memory of the
255  // columns that were used but not the memory of this class member (this should
256  // be fine).
257  void Clear();
258 
259  private:
260  // mutable_column(col) is stored in columns_[mapping_[col]].
261  // The columns_ that can be reused have their index stored in free_columns_.
262  const SparseColumn empty_column_;
264  std::vector<int> free_columns_;
265  std::vector<SparseColumn> columns_;
266  DISALLOW_COPY_AND_ASSIGN(SparseMatrixWithReusableColumnMemory);
267 };
268 
269 // The class that computes either the actual L.U decomposition, or the
270 // permutation P and Q such that P.B.Q^{-1} will have a sparse L.U
271 // decomposition.
272 class Markowitz {
273  public:
275 
276  // Computes the full factorization with P, Q, L and U.
277  //
278  // If the matrix is singular, the returned status will indicate it and the
279  // permutation (col_perm) will contain a maximum non-singular set of columns
280  // of the matrix. Moreover, by adding singleton columns with a one at the rows
281  // such that 'row_perm[row] == kInvalidRow', then the matrix will be
282  // non-singular.
283  ABSL_MUST_USE_RESULT Status
284  ComputeLU(const CompactSparseMatrixView& basis_matrix,
285  RowPermutation* row_perm, ColumnPermutation* col_perm,
286  TriangularMatrix* lower, TriangularMatrix* upper);
287 
288  // Only computes P and Q^{-1}, L and U can be computed later from these
289  // permutations using another algorithm (for instance left-looking L.U). This
290  // may be faster than computing the full L and U at the same time but the
291  // current implementation is not optimized for this.
292  //
293  // It behaves the same as ComputeLU() for singular matrices.
294  //
295  // This function also works with a non-square matrix. It will return a set of
296  // independent columns of maximum size. If all the given columns are
297  // independent, the returned Status will be OK.
298  ABSL_MUST_USE_RESULT Status ComputeRowAndColumnPermutation(
299  const CompactSparseMatrixView& basis_matrix, RowPermutation* row_perm,
300  ColumnPermutation* col_perm);
301 
302  // Releases the memory used by this class.
303  void Clear();
304 
305  // Returns a string containing the statistics for this class.
306  std::string StatString() const { return stats_.StatString(); }
307 
308  // Sets the current parameters.
309  void SetParameters(const GlopParameters& parameters) {
310  parameters_ = parameters;
311  }
312 
313  private:
314  // Statistics about this class.
315  struct Stats : public StatsGroup {
316  Stats()
317  : StatsGroup("Markowitz"),
318  basis_singleton_column_ratio("basis_singleton_column_ratio", this),
319  basis_residual_singleton_column_ratio(
320  "basis_residual_singleton_column_ratio", this),
321  pivots_without_fill_in_ratio("pivots_without_fill_in_ratio", this),
322  degree_two_pivot_columns("degree_two_pivot_columns", this) {}
323  RatioDistribution basis_singleton_column_ratio;
324  RatioDistribution basis_residual_singleton_column_ratio;
325  RatioDistribution pivots_without_fill_in_ratio;
326  RatioDistribution degree_two_pivot_columns;
327  };
328  Stats stats_;
329 
330  // Fast track for singleton columns of the matrix. Fills a part of the row and
331  // column permutation that move these columns in order to form an identity
332  // sub-matrix on the upper left.
333  //
334  // Note(user): Linear programming bases usually have a resonable percentage of
335  // slack columns in them, so this gives a big speedup.
336  void ExtractSingletonColumns(const CompactSparseMatrixView& basis_matrix,
337  RowPermutation* row_perm,
338  ColumnPermutation* col_perm, int* index);
339 
340  // Fast track for columns that form a triangular matrix. This does not find
341  // all of them, but because the column are ordered in the same way they were
342  // ordered at the end of the previous factorization, this is likely to find
343  // quite a few.
344  //
345  // The main gain here is that it avoids taking these columns into account in
346  // InitializeResidualMatrix() and later in RemoveRowFromResidualMatrix().
347  void ExtractResidualSingletonColumns(
348  const CompactSparseMatrixView& basis_matrix, RowPermutation* row_perm,
349  ColumnPermutation* col_perm, int* index);
350 
351  // Helper function for determining if a column is a residual singleton column.
352  // If it is, RowIndex* row contains the index of the single residual edge.
353  bool IsResidualSingletonColumn(const ColumnView& column,
354  const RowPermutation& row_perm, RowIndex* row);
355 
356  // Returns the column of the current residual matrix with an index 'col' in
357  // the initial matrix. We compute it by solving a linear system with the
358  // current lower_ and the last computed column 'col' of a previous residual
359  // matrix. This uses the same algorithm as a left-looking factorization (see
360  // lu_factorization.h for more details).
361  const SparseColumn& ComputeColumn(const RowPermutation& row_perm,
362  ColIndex col);
363 
364  // Finds an entry in the residual matrix with a low Markowitz score and a high
365  // enough magnitude. Returns its Markowitz score and updates the given
366  // pointers.
367  //
368  // We use the strategy of Zlatev, "On some pivotal strategies in Gaussian
369  // elimination by sparse technique" (1980). SIAM J. Numer. Anal. 17 18-30. It
370  // consists of looking for the best pivot in only a few columns (usually 3
371  // or 4) amongst the ones which have the lowest number of entries.
372  //
373  // Amongst the pivots with a minimum Markowitz number, we choose the one
374  // with highest magnitude. This doesn't apply to pivots with a 0 Markowitz
375  // number because all such pivots will have to be taken at some point anyway.
376  int64 FindPivot(const RowPermutation& row_perm,
377  const ColumnPermutation& col_perm, RowIndex* pivot_row,
378  ColIndex* pivot_col, Fractional* pivot_coefficient);
379 
380  // Updates the degree of a given column in the internal structure of the
381  // class.
382  void UpdateDegree(ColIndex col, int degree);
383 
384  // Removes all the coefficients in the residual matrix that are on the given
385  // row or column. In both cases, the pivot row or column is ignored.
386  void RemoveRowFromResidualMatrix(RowIndex pivot_row, ColIndex pivot_col);
387  void RemoveColumnFromResidualMatrix(RowIndex pivot_row, ColIndex pivot_col);
388 
389  // Updates the residual matrix given the pivot position. This is needed if the
390  // pivot row and pivot column both have more than one entry. Otherwise, the
391  // residual matrix can be updated more efficiently by calling one of the
392  // Remove...() functions above.
393  void UpdateResidualMatrix(RowIndex pivot_row, ColIndex pivot_col);
394 
395  // Pointer to the matrix to factorize.
396  CompactSparseMatrixView const* basis_matrix_;
397 
398  // These matrices are transformed during the algorithm into the final L and U
399  // matrices modulo some row and column permutations. Note that the columns of
400  // these matrices stay in the initial order.
401  SparseMatrixWithReusableColumnMemory permuted_lower_;
402  SparseMatrixWithReusableColumnMemory permuted_upper_;
403 
404  // These matrices will hold the final L and U. The are created columns by
405  // columns from left to right, and at the end, their rows are permuted by
406  // ComputeLU() to become triangular.
407  TriangularMatrix lower_;
408  TriangularMatrix upper_;
409 
410  // The columns of permuted_lower_ for which we do need a call to
411  // PermutedLowerSparseSolve(). This speeds up ComputeColumn().
412  DenseBooleanRow permuted_lower_column_needs_solve_;
413 
414  // Contains the non-zero positions of the current residual matrix (the
415  // lower-right square matrix that gets smaller by one row and column at each
416  // Gaussian elimination step).
417  MatrixNonZeroPattern residual_matrix_non_zero_;
418 
419  // Data structure to access the columns by increasing degree.
420  ColumnPriorityQueue col_by_degree_;
421 
422  // True as long as only singleton columns of the residual matrix are used.
423  bool contains_only_singleton_columns_;
424 
425  // Boolean used to know when col_by_degree_ become useful.
426  bool is_col_by_degree_initialized_;
427 
428  // FindPivot() needs to look at the first entries of col_by_degree_, it
429  // temporary put them here before pushing them back to col_by_degree_.
430  std::vector<ColIndex> examined_col_;
431 
432  // Singleton column indices are kept here rather than in col_by_degree_ to
433  // optimize the algorithm: as long as this or singleton_row_ are not empty,
434  // col_by_degree_ do not need to be initialized nor updated.
435  std::vector<ColIndex> singleton_column_;
436 
437  // List of singleton row indices.
438  std::vector<RowIndex> singleton_row_;
439 
440  // Proto holding all the parameters of this algorithm.
441  GlopParameters parameters_;
442 
443  DISALLOW_COPY_AND_ASSIGN(Markowitz);
444 };
445 
446 } // namespace glop
447 } // namespace operations_research
448 
449 #endif // OR_TOOLS_GLOP_MARKOWITZ_H_
#define DCHECK(condition)
Definition: base/logging.h:884
void Reset(int32 max_degree, ColIndex num_cols)
Definition: markowitz.cc:799
void PushOrAdjust(ColIndex col, int32 degree)
Definition: markowitz.cc:807
ABSL_MUST_USE_RESULT Status ComputeLU(const CompactSparseMatrixView &basis_matrix, RowPermutation *row_perm, ColumnPermutation *col_perm, TriangularMatrix *lower, TriangularMatrix *upper)
Definition: markowitz.cc:143
void SetParameters(const GlopParameters &parameters)
Definition: markowitz.h:309
std::string StatString() const
Definition: markowitz.h:306
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
const SparseColumn & column(ColIndex col) const
Definition: markowitz.cc:848
SatParameters parameters
int int32
int64_t int64
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
Permutation< ColIndex > ColumnPermutation
StrictITIVector< ColIndex, bool > DenseBooleanRow
Definition: lp_types.h:302
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