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"
81#include "ortools/glop/parameters.pb.h"
82#include "ortools/glop/status.h"
87#include "ortools/util/stats.h"
88
89namespace operations_research {
90namespace 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.
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().
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.
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