OR-Tools  8.2
primal_edge_norms.cc
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
15
16#include "ortools/base/timer.h"
18
19namespace operations_research {
20namespace glop {
21
23 const VariablesInfo& variables_info,
24 const BasisFactorization& basis_factorization)
25 : compact_matrix_(compact_matrix),
26 variables_info_(variables_info),
27 basis_factorization_(basis_factorization),
28 stats_(),
29 recompute_edge_squared_norms_(true),
30 reset_devex_weights_(true),
31 edge_squared_norms_(),
32 matrix_column_norms_(),
33 devex_weights_(),
34 direction_left_inverse_(),
35 num_operations_(0) {}
36
38 SCOPED_TIME_STAT(&stats_);
39 recompute_edge_squared_norms_ = true;
40 reset_devex_weights_ = true;
41}
42
44 return recompute_edge_squared_norms_;
45}
46
48 if (recompute_edge_squared_norms_) ComputeEdgeSquaredNorms();
49 return edge_squared_norms_;
50}
51
53 if (reset_devex_weights_) ResetDevexWeights();
54 return devex_weights_;
55}
56
58 if (matrix_column_norms_.empty()) ComputeMatrixColumnNorms();
59 return matrix_column_norms_;
60}
61
63 ColIndex entering_col, const ScatteredColumn& direction) {
64 if (!recompute_edge_squared_norms_) {
65 SCOPED_TIME_STAT(&stats_);
66 // Recompute the squared norm of the edge used during this
67 // iteration, i.e. the entering edge. Note the PreciseSquaredNorm()
68 // since it is a small price to pay for an increased precision.
69 const Fractional old_squared_norm = edge_squared_norms_[entering_col];
70 const Fractional precise_squared_norm = 1.0 + PreciseSquaredNorm(direction);
71 edge_squared_norms_[entering_col] = precise_squared_norm;
72
73 const Fractional precise_norm = sqrt(precise_squared_norm);
74 const Fractional estimated_edges_norm_accuracy =
75 (precise_norm - sqrt(old_squared_norm)) / precise_norm;
76 stats_.edges_norm_accuracy.Add(estimated_edges_norm_accuracy);
77 if (std::abs(estimated_edges_norm_accuracy) >
78 parameters_.recompute_edges_norm_threshold()) {
79 VLOG(1) << "Recomputing edge norms: " << sqrt(precise_squared_norm)
80 << " vs " << sqrt(old_squared_norm);
81 recompute_edge_squared_norms_ = true;
82 }
83 }
84}
85
86void PrimalEdgeNorms::UpdateBeforeBasisPivot(ColIndex entering_col,
87 ColIndex leaving_col,
88 RowIndex leaving_row,
89 const ScatteredColumn& direction,
90 UpdateRow* update_row) {
91 SCOPED_TIME_STAT(&stats_);
92 DCHECK_NE(entering_col, leaving_col);
93 if (!recompute_edge_squared_norms_) {
94 update_row->ComputeUpdateRow(leaving_row);
95 ComputeDirectionLeftInverse(entering_col, direction);
96 UpdateEdgeSquaredNorms(entering_col, leaving_col, leaving_row,
97 direction.values, *update_row);
98 }
99 if (!reset_devex_weights_) {
100 // Resets devex weights once in a while. If so, no need to update them
101 // before.
102 ++num_devex_updates_since_reset_;
103 if (num_devex_updates_since_reset_ >
104 parameters_.devex_weights_reset_period()) {
105 reset_devex_weights_ = true;
106 } else {
107 update_row->ComputeUpdateRow(leaving_row);
108 UpdateDevexWeights(entering_col, leaving_col, leaving_row,
109 direction.values, *update_row);
110 }
111 }
112}
113
114void PrimalEdgeNorms::ComputeMatrixColumnNorms() {
115 SCOPED_TIME_STAT(&stats_);
116 matrix_column_norms_.resize(compact_matrix_.num_cols(), 0.0);
117 for (ColIndex col(0); col < compact_matrix_.num_cols(); ++col) {
118 matrix_column_norms_[col] = sqrt(SquaredNorm(compact_matrix_.column(col)));
119 num_operations_ += compact_matrix_.column(col).num_entries().value();
120 }
121}
122
123void PrimalEdgeNorms::ComputeEdgeSquaredNorms() {
124 SCOPED_TIME_STAT(&stats_);
125
126 // Since we will do a lot of inversions, it is better to be as efficient and
127 // precise as possible by refactorizing the basis.
128 DCHECK(basis_factorization_.IsRefactorized());
129 edge_squared_norms_.resize(compact_matrix_.num_cols(), 0.0);
130 for (const ColIndex col : variables_info_.GetIsRelevantBitRow()) {
131 // Note the +1.0 in the squared norm for the component of the edge on the
132 // 'entering_col'.
133 edge_squared_norms_[col] = 1.0 + basis_factorization_.RightSolveSquaredNorm(
134 compact_matrix_.column(col));
135 }
136 recompute_edge_squared_norms_ = false;
137}
138
139// TODO(user): It should be possible to reorganize the code and call this when
140// the value of direction is no longer needed. This will simplify the code and
141// avoid a copy here.
142void PrimalEdgeNorms::ComputeDirectionLeftInverse(
143 ColIndex entering_col, const ScatteredColumn& direction) {
144 SCOPED_TIME_STAT(&stats_);
145
146 // Initialize direction_left_inverse_ to direction. Note the special case when
147 // the non-zero vector is empty which means we don't know and need to use the
148 // dense version.
149 const ColIndex size = RowToColIndex(direction.values.size());
150 const double kThreshold = 0.05 * size.value();
151 if (!direction_left_inverse_.non_zeros.empty() &&
152 (direction_left_inverse_.non_zeros.size() + direction.non_zeros.size() <
153 2 * kThreshold)) {
154 ClearAndResizeVectorWithNonZeros(size, &direction_left_inverse_);
155 for (const auto e : direction) {
156 direction_left_inverse_[RowToColIndex(e.row())] = e.coefficient();
157 }
158 } else {
159 direction_left_inverse_.values = Transpose(direction.values);
160 direction_left_inverse_.non_zeros.clear();
161 }
162
163 if (direction.non_zeros.size() < kThreshold) {
164 direction_left_inverse_.non_zeros = TransposedView(direction).non_zeros;
165 }
166 basis_factorization_.LeftSolve(&direction_left_inverse_);
167
168 // TODO(user): Refactorize if estimated accuracy above a threshold.
169 IF_STATS_ENABLED(stats_.direction_left_inverse_accuracy.Add(
170 compact_matrix_.ColumnScalarProduct(entering_col,
171 direction_left_inverse_.values) -
172 SquaredNorm(direction.values)));
173 IF_STATS_ENABLED(stats_.direction_left_inverse_density.Add(
174 Density(direction_left_inverse_.values)));
175}
176
177// Let new_edge denote the edge of 'col' in the new basis. We want:
178// reduced_costs_[col] = ScalarProduct(new_edge, basic_objective_);
179// edge_squared_norms_[col] = SquaredNorm(new_edge);
180//
181// In order to compute this, we use the formulas:
182// new_leaving_edge = old_entering_edge / divisor.
183// new_edge = old_edge + update_coeff * new_leaving_edge.
184void PrimalEdgeNorms::UpdateEdgeSquaredNorms(ColIndex entering_col,
185 ColIndex leaving_col,
186 RowIndex leaving_row,
187 const DenseColumn& direction,
188 const UpdateRow& update_row) {
189 SCOPED_TIME_STAT(&stats_);
190
191 // 'pivot' is the value of the entering_edge at 'leaving_row'.
192 // The edge of the 'leaving_col' in the new basis is equal to
193 // entering_edge / 'pivot'.
194 const Fractional pivot = -direction[leaving_row];
195 DCHECK_NE(pivot, 0.0);
196
197 // Note that this should be precise because of the call to
198 // TestEnteringEdgeNormPrecision().
199 const Fractional entering_squared_norm = edge_squared_norms_[entering_col];
200 const Fractional leaving_squared_norm =
201 std::max(1.0, entering_squared_norm / Square(pivot));
202
203 int stat_lower_bounded_norms = 0;
204 const Fractional factor = 2.0 / pivot;
205 for (const ColIndex col : update_row.GetNonZeroPositions()) {
206 const Fractional coeff = update_row.GetCoefficient(col);
207 const Fractional scalar_product = compact_matrix_.ColumnScalarProduct(
208 col, direction_left_inverse_.values);
209 num_operations_ += compact_matrix_.column(col).num_entries().value();
210
211 // Update the edge squared norm of this column. Note that the update
212 // formula used is important to maximize the precision. See an explanation
213 // in the dual context in Koberstein's PhD thesis, section 8.2.2.1.
214 edge_squared_norms_[col] +=
215 coeff * (coeff * leaving_squared_norm + factor * scalar_product);
216
217 // Make sure it doesn't go under a known lower bound (TODO(user): ref?).
218 // This way norms are always >= 1.0 .
219 // TODO(user): precompute 1 / Square(pivot) or 1 / pivot? it will be
220 // slightly faster, but may introduce numerical issues. More generally,
221 // this test is only needed in a few cases, so is it worth it?
222 const Fractional lower_bound = 1.0 + Square(coeff / pivot);
223 if (edge_squared_norms_[col] < lower_bound) {
224 edge_squared_norms_[col] = lower_bound;
225 ++stat_lower_bounded_norms;
226 }
227 }
228 edge_squared_norms_[leaving_col] = leaving_squared_norm;
229 stats_.lower_bounded_norms.Add(stat_lower_bounded_norms);
230}
231
232void PrimalEdgeNorms::UpdateDevexWeights(
233 ColIndex entering_col /* index q in the paper */,
234 ColIndex leaving_col /* index p in the paper */, RowIndex leaving_row,
235 const DenseColumn& direction, const UpdateRow& update_row) {
236 SCOPED_TIME_STAT(&stats_);
237
238 // Compared to steepest edge update, the DEVEX weight uses the largest of the
239 // norms of two vectors to approximate the norm of the sum.
240 const Fractional entering_norm = sqrt(PreciseSquaredNorm(direction));
241 const Fractional pivot_magnitude = std::abs(direction[leaving_row]);
242 const Fractional leaving_norm =
243 std::max(1.0, entering_norm / pivot_magnitude);
244 for (const ColIndex col : update_row.GetNonZeroPositions()) {
245 const Fractional coeff = update_row.GetCoefficient(col);
246 const Fractional update_vector_norm = std::abs(coeff) * leaving_norm;
247 devex_weights_[col] = std::max(devex_weights_[col], update_vector_norm);
248 }
249 devex_weights_[leaving_col] = leaving_norm;
250}
251
252void PrimalEdgeNorms::ResetDevexWeights() {
253 SCOPED_TIME_STAT(&stats_);
254 if (parameters_.initialize_devex_with_column_norms()) {
255 devex_weights_ = GetMatrixColumnNorms();
256 } else {
257 devex_weights_.assign(compact_matrix_.num_cols(), 1.0);
258 }
259 num_devex_updates_since_reset_ = 0;
260 reset_devex_weights_ = false;
261}
262
263} // namespace glop
264} // namespace operations_research
int64 max
Definition: alldiff_cst.cc:139
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define DCHECK(condition)
Definition: base/logging.h:884
#define VLOG(verboselevel)
Definition: base/logging.h:978
bool empty() const
Fractional RightSolveSquaredNorm(const ColumnView &a) const
Fractional ColumnScalarProduct(ColIndex col, const DenseRow &vector) const
Definition: sparse.h:382
ColumnView column(ColIndex col) const
Definition: sparse.h:364
PrimalEdgeNorms(const CompactSparseMatrix &compact_matrix, const VariablesInfo &variables_info, const BasisFactorization &basis_factorization)
void UpdateBeforeBasisPivot(ColIndex entering_col, ColIndex leaving_col, RowIndex leaving_row, const ScatteredColumn &direction, UpdateRow *update_row)
void TestEnteringEdgeNormPrecision(ColIndex entering_col, const ScatteredColumn &direction)
void assign(IntType size, const T &v)
Definition: lp_types.h:274
void ComputeUpdateRow(RowIndex leaving_row)
Definition: update_row.cc:78
const DenseBitRow & GetIsRelevantBitRow() const
ColIndex col
Definition: markowitz.cc:176
Fractional PreciseSquaredNorm(const SparseColumn &v)
Fractional SquaredNorm(const SparseColumn &v)
Fractional Square(Fractional f)
double Density(const DenseRow &row)
ColIndex RowToColIndex(RowIndex row)
Definition: lp_types.h:48
const DenseRow & Transpose(const DenseColumn &col)
const ScatteredRow & TransposedView(const ScatteredColumn &c)
void ClearAndResizeVectorWithNonZeros(IndexType size, ScatteredRowOrCol *v)
StrictITIVector< RowIndex, Fractional > DenseColumn
Definition: lp_types.h:328
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
#define IF_STATS_ENABLED(instructions)
Definition: stats.h:437
#define SCOPED_TIME_STAT(stats)
Definition: stats.h:438
StrictITIVector< Index, Fractional > values