OR-Tools  8.2
lu_factorization.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 #ifndef OR_TOOLS_GLOP_LU_FACTORIZATION_H_
15 #define OR_TOOLS_GLOP_LU_FACTORIZATION_H_
16 
17 #include "ortools/glop/markowitz.h"
18 #include "ortools/glop/parameters.pb.h"
19 #include "ortools/glop/status.h"
23 #include "ortools/lp_data/sparse.h"
25 #include "ortools/util/stats.h"
26 
27 namespace operations_research {
28 namespace glop {
29 
30 // An LU-Factorization class encapsulating the LU factorization data and
31 // algorithms. The actual algorithm is in markowitz.h and .cc. This class holds
32 // all the Solve() functions that deal with the permutations and the L and U
33 // factors once they are computed.
35  public:
37 
38  // Returns true if the LuFactorization is a factorization of the identity
39  // matrix. In this state, all the Solve() functions will work for any
40  // vector dimension.
41  bool IsIdentityFactorization() { return is_identity_factorization_; }
42 
43  // Clears internal data structure and reset this class to the factorization
44  // of an identity matrix.
45  void Clear();
46 
47  // Computes an LU-decomposition for a given matrix B. If for some reason,
48  // there was an error, then the factorization is reset to the one of the
49  // identity matrix, and an error is reported.
50  //
51  // Note(user): Since a client must use the result, there is little chance of
52  // it being confused by this revert to identity factorization behavior. The
53  // reason behind it is that this way, calling any public function of this
54  // class will never cause a crash of the program.
55  ABSL_MUST_USE_RESULT Status
56  ComputeFactorization(const CompactSparseMatrixView& compact_matrix);
57 
58  // Returns the column permutation used by the LU factorization.
59  const ColumnPermutation& GetColumnPermutation() const { return col_perm_; }
60 
61  // Sets the column permutation to the identity permutation. The idea is that
62  // the column permutation can be incorporated in the basis RowToColMapping,
63  // and once this is done, then a client can call this and effectively remove
64  // the need for a column permutation on each solve.
66  col_perm_.clear();
67  inverse_col_perm_.clear();
68  }
69 
70  // Solves 'B.x = b', x initially contains b, and is replaced by 'B^{-1}.b'.
71  // Since P.B.Q^{-1} = L.U, we have B = P^{-1}.L.U.Q.
72  // 1/ Solve P^{-1}.y = b for y by computing y = P.b,
73  // 2/ solve L.z = y for z,
74  // 3/ solve U.t = z for t,
75  // 4/ finally solve Q.x = t, by computing x = Q^{-1}.t.
76  void RightSolve(DenseColumn* x) const;
77 
78  // Solves 'y.B = r', y initially contains r, and is replaced by r.B^{-1}.
79  // Internally, it takes x = y^T, b = r^T and solves B^T.x = b.
80  // We have P.B.Q^{-1} = P.B.Q^T = L.U, thus (L.U)^T = Q.B^T.P^T.
81  // Therefore B^T = Q^{-1}.U^T.L^T.P^T.P^{-1} = Q^{-1}.U^T.L^T.P
82  // The procedure is thus:
83  // 1/ Solve Q^{-1}.y = b for y, by computing y = Q.b,
84  // 2/ solve U^T.z = y for z,
85  // 3/ solve L^T.t = z for t,
86  // 4/ finally, solve P.x = t for x by computing x = P^{-1}.t.
87  void LeftSolve(DenseRow* y) const;
88 
89  // More fine-grained right/left solve functions that may exploit the initial
90  // non-zeros of the input vector if non-empty. Note that a solve involving L
91  // actually solves P^{-1}.L and a solve involving U actually solves U.Q. To
92  // solve a system with the initial matrix B, one needs to call:
93  // - RightSolveL() and then RightSolveU() for a right solve (B.x = initial x).
94  // - LeftSolveU() and then LeftSolveL() for a left solve (y.B = initial y).
97  void LeftSolveUWithNonZeros(ScatteredRow* y) const;
98 
99  // Specialized version of LeftSolveL() that may exploit the initial non_zeros
100  // of y if it is non empty. Moreover, if result_before_permutation is not
101  // NULL, it might be filled with the result just before row_perm_ is applied
102  // to it and true is returned. If result_before_permutation is not filled,
103  // then false is returned.
105  ScatteredColumn* result_before_permutation) const;
106  void LeftSolveLWithNonZeros(ScatteredRow* y) const;
107 
108  // Specialized version of RightSolveLWithNonZeros() that takes a SparseColumn
109  // or a ScatteredColumn as input. non_zeros will either be cleared or set to
110  // the non zeros of the result. Important: the output x must be of the correct
111  // size and all zero.
112  void RightSolveLForColumnView(const ColumnView& b, ScatteredColumn* x) const;
114  ScatteredColumn* x) const;
115 
116  // Specialized version of RightSolveLWithNonZeros() where x is originaly equal
117  // to 'a' permuted by row_perm_. Note that 'a' is only used for DCHECK.
119  ScatteredColumn* x) const;
120 
121  // Specialized version of LeftSolveU() for an unit right-hand side.
122  // non_zeros will either be cleared or set to the non zeros of the results.
123  // It also returns the value of col permuted by Q (which is the position
124  // of the unit-vector rhs in the solve system: y.U = rhs).
125  // Important: the output y must be of the correct size and all zero.
126  ColIndex LeftSolveUForUnitRow(ColIndex col, ScatteredRow* y) const;
127 
128  // Returns the given column of U.
129  // It will only be valid until the next call to GetColumnOfU().
130  const SparseColumn& GetColumnOfU(ColIndex col) const;
131 
132  // Returns the norm of B^{-1}.a
134 
135  // Returns the norm of (B^T)^{-1}.e_row where e is an unit vector.
136  Fractional DualEdgeSquaredNorm(RowIndex row) const;
137 
138  // The fill-in of the LU-factorization is defined as the sum of the number
139  // of entries of both the lower- and upper-triangular matrices L and U minus
140  // the number of entries in the initial matrix B.
141  //
142  // This returns the number of entries in lower + upper as the percentage of
143  // the number of entries in B.
144  double GetFillInPercentage(const CompactSparseMatrixView& matrix) const;
145 
146  // Returns the number of entries in L + U.
147  // If the factorization is the identity, this returns 0.
148  EntryIndex NumberOfEntries() const;
149 
150  // Computes the determinant of the input matrix B.
151  // Since P.B.Q^{-1} = L.U, det(P) * det(B) * det(Q^{-1}) = det(L) * det(U).
152  // det(L) = 1 since L is a lower-triangular matrix with 1 on the diagonal.
153  // det(P) = +1 or -1 (by definition it is the sign of the permutation P)
154  // det(Q^{-1}) = +1 or -1 (the sign of the permutation Q^{-1})
155  // Finally det(U) = product of the diagonal elements of U, since U is an
156  // upper-triangular matrix.
157  // Taking all this into account:
158  // det(B) = sign(P) * sign(Q^{-1}) * prod_i u_ii .
160 
161  // Computes the 1-norm of the inverse of the input matrix B.
162  // For this we iteratively solve B.x = e_j, where e_j is the jth unit vector.
163  // The result of this computation is the jth column of B^-1.
164  // The 1-norm |B| is defined as max_j sum_i |a_ij|
165  // http://en.wikipedia.org/wiki/Matrix_norm
167 
168  // Computes the infinity-norm of the inverse of the input matrix B.
169  // The infinity-norm |B| is defined as max_i sum_j |a_ij|
170  // http://en.wikipedia.org/wiki/Matrix_norm
172 
173  // Computes the condition number of the input matrix B.
174  // For a given norm, this is the matrix norm times the norm of its inverse.
175  //
176  // Note that because the LuFactorization class does not keep the
177  // non-factorized matrix in memory, it needs to be passed to these functions.
178  // It is up to the client to pass exactly the same matrix as the one used
179  // for ComputeFactorization().
180  //
181  // TODO(user): separate this from LuFactorization.
183  const CompactSparseMatrixView& matrix) const;
185  const CompactSparseMatrixView& matrix) const;
187 
188  // Sets the current parameters.
189  void SetParameters(const GlopParameters& parameters) {
190  parameters_ = parameters;
191  markowitz_.SetParameters(parameters);
192  }
193 
194  // Returns a string containing the statistics for this class.
195  std::string StatString() const {
196  return stats_.StatString() + markowitz_.StatString();
197  }
198 
199  // This is only used for testing and in debug mode.
200  // TODO(user): avoid the matrix conversion by multiplying TriangularMatrix
201  // directly.
202  void ComputeLowerTimesUpper(SparseMatrix* product) const {
203  SparseMatrix temp_lower, temp_upper;
204  lower_.CopyToSparseMatrix(&temp_lower);
205  upper_.CopyToSparseMatrix(&temp_upper);
206  product->PopulateFromProduct(temp_lower, temp_upper);
207  }
208 
209  // Visible for testing.
210  const RowPermutation& row_perm() const { return row_perm_; }
212  return inverse_col_perm_;
213  }
214 
215  private:
216  // Statistics about this class.
217  struct Stats : public StatsGroup {
218  Stats()
219  : StatsGroup("LuFactorization"),
220  basis_num_entries("basis_num_entries", this),
221  lu_fill_in("lu_fill_in", this) {}
222  IntegerDistribution basis_num_entries;
223  RatioDistribution lu_fill_in;
224  };
225 
226  // Internal function used in the left solve functions.
227  void LeftSolveScratchpad() const;
228 
229  // Internal function used in the right solve functions
230  template <typename Column>
231  void RightSolveLInternal(const Column& b, ScatteredColumn* x) const;
232 
233  // Fills transpose_upper_ from upper_.
234  void ComputeTransposeUpper();
235 
236  // transpose_lower_ is only needed when we compute dual norms.
237  void ComputeTransposeLower() const;
238 
239  // Computes R = P.B.Q^{-1} - L.U and returns false if the largest magnitude of
240  // the coefficients of P.B.Q^{-1} - L.U is greater than tolerance.
241  bool CheckFactorization(const CompactSparseMatrixView& matrix,
242  Fractional tolerance) const;
243 
244  // Special case where we have nothing to do. This happens at the beginning
245  // when we start the problem with an all-slack basis and gives a good speedup
246  // on really easy problems. It is initially true and set to true each time we
247  // call Clear(). We set it to false if a call to ComputeFactorization()
248  // succeeds.
249  bool is_identity_factorization_;
250 
251  // The triangular factor L and U (and its transpose).
252  TriangularMatrix lower_;
253  TriangularMatrix upper_;
254  TriangularMatrix transpose_upper_;
255 
256  // The transpose of lower_. It is just used by DualEdgeSquaredNorm()
257  // and mutable so it can be lazily initialized.
258  mutable TriangularMatrix transpose_lower_;
259 
260  // The column permutation Q and its inverse Q^{-1} in P.B.Q^{-1} = L.U.
261  ColumnPermutation col_perm_;
262  ColumnPermutation inverse_col_perm_;
263 
264  // The row permutation P and its inverse P^{-1} in P.B.Q^{-1} = L.U.
265  RowPermutation row_perm_;
266  RowPermutation inverse_row_perm_;
267 
268  // Temporary storage used by LeftSolve()/RightSolve().
269  mutable DenseColumn dense_column_scratchpad_;
270 
271  // Temporary storage used by GetColumnOfU().
272  mutable SparseColumn column_of_upper_;
273 
274  // Same as dense_column_scratchpad_ but this vector is always reset to zero by
275  // the functions that use it. non_zero_rows_ is used to track the
276  // non_zero_rows_ position of dense_column_scratchpad_.
277  mutable DenseColumn dense_zero_scratchpad_;
278  mutable std::vector<RowIndex> non_zero_rows_;
279 
280  // Statistics, mutable so const functions can still update it.
281  mutable Stats stats_;
282 
283  // Proto holding all the parameters of this algorithm.
284  GlopParameters parameters_;
285 
286  // The class doing the Markowitz LU factorization.
287  Markowitz markowitz_;
288 
289  DISALLOW_COPY_AND_ASSIGN(LuFactorization);
290 };
291 
292 } // namespace glop
293 } // namespace operations_research
294 #endif // OR_TOOLS_GLOP_LU_FACTORIZATION_H_
void LeftSolveUWithNonZeros(ScatteredRow *y) const
const SparseColumn & GetColumnOfU(ColIndex col) const
void RightSolveLForColumnView(const ColumnView &b, ScatteredColumn *x) const
void RightSolveLWithPermutedInput(const DenseColumn &a, ScatteredColumn *x) const
double GetFillInPercentage(const CompactSparseMatrixView &matrix) const
Fractional RightSolveSquaredNorm(const ColumnView &a) const
const ColumnPermutation & GetColumnPermutation() const
void RightSolveUWithNonZeros(ScatteredColumn *x) const
const RowPermutation & row_perm() const
bool LeftSolveLWithNonZeros(ScatteredRow *y, ScatteredColumn *result_before_permutation) const
const ColumnPermutation & inverse_col_perm() const
ColIndex LeftSolveUForUnitRow(ColIndex col, ScatteredRow *y) const
Fractional DualEdgeSquaredNorm(RowIndex row) const
void RightSolveLForScatteredColumn(const ScatteredColumn &b, ScatteredColumn *x) const
void RightSolveLWithNonZeros(ScatteredColumn *x) const
Fractional ComputeInfinityNormConditionNumber(const CompactSparseMatrixView &matrix) const
void ComputeLowerTimesUpper(SparseMatrix *product) const
ABSL_MUST_USE_RESULT Status ComputeFactorization(const CompactSparseMatrixView &compact_matrix)
void SetParameters(const GlopParameters &parameters)
Fractional ComputeOneNormConditionNumber(const CompactSparseMatrixView &matrix) const
void SetParameters(const GlopParameters &parameters)
Definition: markowitz.h:309
std::string StatString() const
Definition: markowitz.h:306
void PopulateFromProduct(const SparseMatrix &a, const SparseMatrix &b)
Definition: sparse.cc:250
void CopyToSparseMatrix(SparseMatrix *output) const
Definition: sparse.cc:730
SatParameters parameters
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
Permutation< ColIndex > ColumnPermutation
StrictITIVector< RowIndex, Fractional > DenseColumn
Definition: lp_types.h:328
Permutation< RowIndex > RowPermutation
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...