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
18#include "ortools/glop/parameters.pb.h"
19#include "ortools/glop/status.h"
25#include "ortools/util/stats.h"
26
27namespace operations_research {
28namespace 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).
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;
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.
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.
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
const ColumnPermutation & inverse_col_perm() const
void RightSolveLWithPermutedInput(const DenseColumn &a, ScatteredColumn *x) const
double GetFillInPercentage(const CompactSparseMatrixView &matrix) const
Fractional RightSolveSquaredNorm(const ColumnView &a) const
void RightSolveUWithNonZeros(ScatteredColumn *x) const
const ColumnPermutation & GetColumnPermutation() const
bool LeftSolveLWithNonZeros(ScatteredRow *y, ScatteredColumn *result_before_permutation) 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
const RowPermutation & row_perm() 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...