OR-Tools  8.2
basis_representation.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_BASIS_REPRESENTATION_H_
15#define OR_TOOLS_GLOP_BASIS_REPRESENTATION_H_
16
19#include "ortools/glop/parameters.pb.h"
21#include "ortools/glop/status.h"
26#include "ortools/util/stats.h"
27
28namespace operations_research {
29namespace glop {
30
31// An eta matrix E corresponds to the identity matrix except for one column e of
32// index j. In particular, B.E is the matrix of the new basis obtained from B by
33// replacing the j-th vector of B by B.e, note that this is exactly what happens
34// during a "pivot" of the current basis in the simplex algorithm.
35//
36// E = [ 1 ... 0 e_0 0 ... 0
37// ... ... ... ... ... ... ...
38// 0 ... 1 e_{j-1} 0 ... 0
39// 0 ... 0 e_j 0 ... 0
40// 0 ... 0 e_{j+1} 1 ... 0
41// ... ... ... ... ... ... ...
42// 0 ... 0 e_{n-1} 0 ... 1 ]
43//
44// The inverse of the eta matrix is:
45// E^{-1} = [ 1 ... 0 -e_0/e_j 0 ... 0
46// ... ... ... ... ... ... ...
47// 0 ... 1 -e_{j-1}/e_j 0 ... 0
48// 0 ... 0 1/e_j 0 ... 0
49// 0 ... 0 -e_{j+1}/e_j 1 ... 0
50// ... ... ... ... ... ... ...
51// 0 ... 0 -e_{n-1}/e_j 0 ... 1 ]
52class EtaMatrix {
53 public:
54 EtaMatrix(ColIndex eta_col, const ScatteredColumn& direction);
55 virtual ~EtaMatrix();
56
57 // Solves the system y.E = c, 'c' beeing the initial value of 'y'.
58 // Then y = c.E^{-1}, so y is equal to c except for
59 // y_j = (c_j - \sum_{i != j}{c_i * e_i}) / e_j.
60 void LeftSolve(DenseRow* y) const;
61
62 // Same as LeftSolve(), but 'pos' contains the non-zero positions of c. The
63 // order of the positions is not important, but there must be no duplicates.
64 // The values not in 'pos' are not used. If eta_col_ was not already in 'pos',
65 // it is added.
66 void SparseLeftSolve(DenseRow* y, ColIndexVector* pos) const;
67
68 // Solves the system E.d = a, 'a' beeing the initial value of 'd'.
69 // Then d = E^{-1}.a = [ a_0 - e_0 * a_j / e_j
70 // ...
71 // a_{j-1} - e_{j-1} * a_j / e_j
72 // a_j / e_j
73 // a_{j+1} - e_{j+1} * a_j / e_j
74 // ...
75 // a_{n-1} - e_{n-1} * a_j / e_j ]
76 void RightSolve(DenseColumn* d) const;
77
78 private:
79 // Internal RightSolve() and LeftSolve() implementations using either the
80 // dense or the sparse representation of the eta vector.
81 void LeftSolveWithDenseEta(DenseRow* y) const;
82 void LeftSolveWithSparseEta(DenseRow* y) const;
83 void RightSolveWithDenseEta(DenseColumn* d) const;
84 void RightSolveWithSparseEta(DenseColumn* d) const;
85
86 // If an eta vector density is smaller than this threshold, we use the
87 // sparse version of the Solve() functions rather than the dense version.
88 // TODO(user): Detect automatically a good parameter? 0.5 is a good value on
89 // the Netlib (I only did a few experiments though). Note that in the future
90 // we may not even keep the dense representation at all.
91 static const Fractional kSparseThreshold;
92
93 const ColIndex eta_col_;
94 const Fractional eta_col_coefficient_;
95
96 // Note that to optimize solves, the position eta_col_ is set to 0.0 and
97 // stored in eta_col_coefficient_ instead.
98 DenseColumn eta_coeff_;
99 SparseColumn sparse_eta_coeff_;
100
101 DISALLOW_COPY_AND_ASSIGN(EtaMatrix);
102};
103
104// An eta factorization corresponds to the product of k eta matrices,
105// i.e. E = E_0.E_1. ... .E_{k-1}
106// It is used to solve two systems:
107// - E.d = a (where a is usually the entering column).
108// - y.E = c (where c is usually the objective row).
110 public:
112 virtual ~EtaFactorization();
113
114 // Deletes all eta matrices.
115 void Clear();
116
117 // Updates the eta factorization, i.e. adds the new eta matrix defined by
118 // the leaving variable and the corresponding eta column.
119 void Update(ColIndex entering_col, RowIndex leaving_variable_row,
120 const ScatteredColumn& direction);
121
122 // Left solves all systems from right to left, i.e. y_i = y_{i+1}.(E_i)^{-1}
123 void LeftSolve(DenseRow* y) const;
124
125 // Same as LeftSolve(), but 'pos' contains the non-zero positions of c. The
126 // order of the positions is not important, but there must be no duplicates.
127 // The values not in 'pos' are not used. If eta_col_ was not already in 'pos',
128 // it is added.
129 void SparseLeftSolve(DenseRow* y, ColIndexVector* pos) const;
130
131 // Right solves all systems from left to right, i.e. E_i.d_{i+1} = d_i
132 void RightSolve(DenseColumn* d) const;
133
134 private:
135 std::vector<EtaMatrix*> eta_matrix_;
136
137 DISALLOW_COPY_AND_ASSIGN(EtaFactorization);
138};
139
140// A basis factorization is the product of an eta factorization and
141// a L.U decomposition, i.e. B = L.U.E_0.E_1. ... .E_{k-1}
142// It is used to solve two systems:
143// - B.d = a where a is the entering column.
144// - y.B = c where c is the objective row.
145//
146// To speed-up and improve stability the factorization is refactorized at least
147// every 'refactorization_period' updates.
148//
149// This class does not take ownership of the underlying matrix and basis, and
150// thus they must outlive this class (and keep the same address in memory).
152 public:
153 BasisFactorization(const CompactSparseMatrix* compact_matrix,
154 const RowToColMapping* basis);
155 virtual ~BasisFactorization();
156
157 // Sets the parameters for this component.
158 void SetParameters(const GlopParameters& parameters) {
159 max_num_updates_ = parameters.basis_refactorization_period();
160 use_middle_product_form_update_ =
161 parameters.use_middle_product_form_update();
162 parameters_ = parameters;
163 lu_factorization_.SetParameters(parameters);
164 }
165
166 // Returns the column permutation used by the LU factorization.
167 // This call only makes sense if the basis was just refactorized.
170 return lu_factorization_.GetColumnPermutation();
171 }
172
173 // Sets the column permutation used by the LU factorization to the identity.
174 // Hense the Solve() results will be computed without this permutation.
175 // This call only makes sense if the basis was just refactorized.
178 lu_factorization_.SetColumnPermutationToIdentity();
179 }
180
181 // Clears the factorization and resets it to an identity matrix of size given
182 // by matrix_.num_rows().
183 void Clear();
184
185 // Clears the factorization and initializes the class using the current
186 // matrix_ and basis_. This is fast if IsIdentityBasis() is true, otherwise
187 // it will trigger a refactorization and will return an error if the matrix
188 // could not be factorized.
189 ABSL_MUST_USE_RESULT Status Initialize();
190
191 // Return the number of rows in the basis.
192 RowIndex GetNumberOfRows() const { return compact_matrix_.num_rows(); }
193
194 // Clears eta factorization and refactorizes LU.
195 // Nothing happens if this is called on an already refactorized basis.
196 // Returns an error if the matrix could not be factorized: i.e. not a basis.
197 ABSL_MUST_USE_RESULT Status Refactorize();
198
199 // Like Refactorize(), but do it even if IsRefactorized() is true.
200 // Call this if the underlying basis_ changed and Update() wasn't called.
201 ABSL_MUST_USE_RESULT Status ForceRefactorization();
202
203 // Returns true if the factorization was just recomputed.
204 bool IsRefactorized() const;
205
206 // Updates the factorization. The 'eta' column will be modified with a swap to
207 // avoid a copy (only if the standard eta update is used). Returns an error if
208 // the matrix could not be factorized: i.e. not a basis.
209 ABSL_MUST_USE_RESULT Status Update(ColIndex entering_col,
210 RowIndex leaving_variable_row,
211 const ScatteredColumn& direction);
212
213 // Left solves the system y.B = rhs, where y initialy contains rhs.
214 void LeftSolve(ScatteredRow* y) const;
215
216 // Left solves the system y.B = e_j, where e_j has only 1 non-zero
217 // coefficient of value 1.0 at position 'j'.
218 void LeftSolveForUnitRow(ColIndex j, ScatteredRow* y) const;
219
220 // Same as LeftSolveForUnitRow() but does not update any internal data.
221 void TemporaryLeftSolveForUnitRow(ColIndex j, ScatteredRow* y) const;
222
223 // Right solves the system B.d = a where the input is the initial value of d.
224 void RightSolve(ScatteredColumn* d) const;
225
226 // Same as RightSolve() for matrix.column(col). This also exploits its
227 // sparsity.
228 void RightSolveForProblemColumn(ColIndex col, ScatteredColumn* d) const;
229
230 // Specialized version for ComputeTau() in DualEdgeNorms. This reuses an
231 // intermediate result of the last LeftSolveForUnitRow() in order to save a
232 // permutation if it is available. Note that the input 'a' should always be
233 // equal to the last result of LeftSolveForUnitRow() and will be used for a
234 // DCHECK() or if the intermediate result wasn't kept.
235 const DenseColumn& RightSolveForTau(const ScatteredColumn& a) const;
236
237 // Returns the norm of B^{-1}.a, this is a specific function because
238 // it is a bit faster and it avoids polluting the stats of RightSolve().
239 // It can be called only when IsRefactorized() is true.
241
242 // Returns the norm of (B^T)^{-1}.e_row where e is an unit vector.
243 // This is a bit faster and avoids polluting the stats of LeftSolve().
244 // It can be called only when IsRefactorized() is true.
245 Fractional DualEdgeSquaredNorm(RowIndex row) const;
246
247 // Computes the condition number of B.
248 // For a given norm, this is the matrix norm times the norm of its inverse.
249 // A condition number greater than 1E7 will lead to precision problems.
253
254 // Computes the 1-norm of B.
255 // The 1-norm |A| is defined as max_j sum_i |a_ij|
256 // http://en.wikipedia.org/wiki/Matrix_norm
258
259 // Computes the infinity-norm of B.
260 // The infinity-norm |A| is defined as max_i sum_j |a_ij|
261 // http://en.wikipedia.org/wiki/Matrix_norm
263
264 // Computes the 1-norm of the inverse of B.
265 // For this we iteratively solve B.x = e_j, where e_j is the jth unit vector.
266 // The result of this computation is the jth column of B^-1.
268
269 // Computes the infinity-norm of the inverse of B.
271
272 // Stats related function.
273 // Note that ResetStats() could be const, but until needed it is not to
274 // prevent anyone holding a const BasisFactorization& to call it.
275 std::string StatString() const {
276 return stats_.StatString() + lu_factorization_.StatString();
277 }
278 void ResetStats() { stats_.Reset(); }
279
280 // The deterministic time used by this class. It is incremented for each
281 // solve and each factorization.
282 double DeterministicTime() const;
283
284 private:
285 // Return true if the submatrix of matrix_ given by basis_ is exactly the
286 // identity (without permutation).
287 bool IsIdentityBasis() const;
288
289 // Updates the factorization using the middle product form update.
290 // Qi Huangfu, J. A. Julian Hall, "Novel update techniques for the revised
291 // simplex method", 28 january 2013, Technical Report ERGO-13-0001
292 ABSL_MUST_USE_RESULT Status
293 MiddleProductFormUpdate(ColIndex entering_col, RowIndex leaving_variable_row);
294
295 // Increases the deterministic time for a solve operation with a vector having
296 // this number of non-zero entries (it can be an approximation).
297 void BumpDeterministicTimeForSolve(int num_entries) const;
298
299 // Stats about this class.
300 struct Stats : public StatsGroup {
301 Stats()
302 : StatsGroup("BasisFactorization"),
303 refactorization_interval("refactorization_interval", this) {}
304 IntegerDistribution refactorization_interval;
305 };
306
307 // Mutable because we track the running time of const method like
308 // RightSolve() and LeftSolve().
309 mutable Stats stats_;
310 GlopParameters parameters_;
311
312 // References to the basis subpart of the linear program matrix.
313 const CompactSparseMatrix& compact_matrix_;
314 const RowToColMapping& basis_;
315
316 // Middle form product update factorization and scratchpad_ used to construct
317 // new rank one matrices.
318 RankOneUpdateFactorization rank_one_factorization_;
319 mutable DenseColumn scratchpad_;
320 mutable std::vector<RowIndex> scratchpad_non_zeros_;
321
322 // This is used by RightSolveForTau(). It holds an intermediate result from
323 // the last LeftSolveForUnitRow() and also the final result of
324 // RightSolveForTau().
325 mutable ScatteredColumn tau_;
326
327 // Booleans controlling the interaction between LeftSolveForUnitRow() that may
328 // or may not keep its intermediate results for the optimized
329 // RightSolveForTau().
330 //
331 // tau_computation_can_be_optimized_ will be true iff LeftSolveForUnitRow()
332 // kept its intermediate result when it was called and the factorization
333 // didn't change since then. If it is true, then RightSolveForTau() can use
334 // this result for a faster computation.
335 //
336 // tau_is_computed_ is used as an heuristic by LeftSolveForUnitRow() to decide
337 // if it is worth keeping its intermediate result (which is sligthly slower).
338 // It is simply set to true by RightSolveForTau() and to false by
339 // LeftSolveForUnitRow(), this way the optimization will automatically switch
340 // itself on when switching from the primal simplex (where RightSolveForTau()
341 // is never called) to the dual where it is called after each
342 // LeftSolveForUnitRow(), and back off again in the other direction.
343 mutable bool tau_computation_can_be_optimized_;
344 mutable bool tau_is_computed_;
345
346 // Data structure to store partial solve results for the middle form product
347 // update. See LeftSolveForUnitRow() and RightSolveForProblemColumn(). We use
348 // two CompactSparseMatrix to have a better cache behavior when solving with
349 // the rank_one_factorization_.
350 mutable CompactSparseMatrix storage_;
351 mutable CompactSparseMatrix right_storage_;
352 mutable ColMapping left_pool_mapping_;
353 mutable ColMapping right_pool_mapping_;
354
355 bool use_middle_product_form_update_;
356 int max_num_updates_;
357 int num_updates_;
358 EtaFactorization eta_factorization_;
359 LuFactorization lu_factorization_;
360
361 // mutable because the Solve() functions are const but need to update this.
362 mutable double deterministic_time_;
363
364 DISALLOW_COPY_AND_ASSIGN(BasisFactorization);
365};
366
367} // namespace glop
368} // namespace operations_research
369
370#endif // OR_TOOLS_GLOP_BASIS_REPRESENTATION_H_
#define DCHECK(condition)
Definition: base/logging.h:884
BasisFactorization(const CompactSparseMatrix *compact_matrix, const RowToColMapping *basis)
const DenseColumn & RightSolveForTau(const ScatteredColumn &a) const
void LeftSolveForUnitRow(ColIndex j, ScatteredRow *y) const
Fractional RightSolveSquaredNorm(const ColumnView &a) const
const ColumnPermutation & GetColumnPermutation() const
void TemporaryLeftSolveForUnitRow(ColIndex j, ScatteredRow *y) const
ABSL_MUST_USE_RESULT Status Update(ColIndex entering_col, RowIndex leaving_variable_row, const ScatteredColumn &direction)
Fractional DualEdgeSquaredNorm(RowIndex row) const
void RightSolveForProblemColumn(ColIndex col, ScatteredColumn *d) const
void SetParameters(const GlopParameters &parameters)
void SparseLeftSolve(DenseRow *y, ColIndexVector *pos) const
void Update(ColIndex entering_col, RowIndex leaving_variable_row, const ScatteredColumn &direction)
EtaMatrix(ColIndex eta_col, const ScatteredColumn &direction)
void SparseLeftSolve(DenseRow *y, ColIndexVector *pos) const
const ColumnPermutation & GetColumnPermutation() const
void SetParameters(const GlopParameters &parameters)
SatParameters parameters
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
std::vector< ColIndex > ColIndexVector
Definition: lp_types.h:308
StrictITIVector< RowIndex, ColIndex > RowToColMapping
Definition: lp_types.h:342
StrictITIVector< RowIndex, Fractional > DenseColumn
Definition: lp_types.h:328
StrictITIVector< ColIndex, ColIndex > ColMapping
Definition: lp_types.h:305
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
EntryIndex num_entries