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 
19 namespace operations_research {
20 namespace 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 
86 void 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 
114 void 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 
123 void 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.
142 void 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.
184 void 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 
232 void 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 
252 void 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 Square(Fractional f)
Fractional PreciseSquaredNorm(const SparseColumn &v)
Fractional SquaredNorm(const SparseColumn &v)
double Density(const DenseRow &row)
ColIndex RowToColIndex(RowIndex row)
Definition: lp_types.h:48
void ClearAndResizeVectorWithNonZeros(IndexType size, ScatteredRowOrCol *v)
const DenseRow & Transpose(const DenseColumn &col)
StrictITIVector< RowIndex, Fractional > DenseColumn
Definition: lp_types.h:328
const ScatteredRow & TransposedView(const ScatteredColumn &c)
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