OR-Tools  8.2
dual_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 
17 
18 namespace operations_research {
19 namespace glop {
20 
22  : basis_factorization_(basis_factorization),
23  recompute_edge_squared_norms_(true) {}
24 
26  return recompute_edge_squared_norms_;
27 }
28 
29 void DualEdgeNorms::Clear() { recompute_edge_squared_norms_ = true; }
30 
31 void DualEdgeNorms::ResizeOnNewRows(RowIndex new_size) {
32  edge_squared_norms_.resize(new_size, 1.0);
33 }
34 
36  if (recompute_edge_squared_norms_) ComputeEdgeSquaredNorms();
37  return edge_squared_norms_;
38 }
39 
41  const ColumnPermutation& col_perm) {
42  if (recompute_edge_squared_norms_) return;
43  ApplyColumnPermutationToRowIndexedVector(col_perm, &edge_squared_norms_);
44 }
45 
47  ColIndex entering_col, RowIndex leaving_row,
48  const ScatteredColumn& direction,
49  const ScatteredRow& unit_row_left_inverse) {
50  // No need to update if we will recompute it from scratch later.
51  if (recompute_edge_squared_norms_) return;
52  const DenseColumn& tau = ComputeTau(TransposedView(unit_row_left_inverse));
53  SCOPED_TIME_STAT(&stats_);
54 
55  // ||unit_row_left_inverse||^2 is the same as
56  // edge_squared_norms_[leaving_row], but with a better precision. If the
57  // difference between the two is too large, we trigger a full recomputation.
58  //
59  // Note that we use the PreciseSquaredNorms() because it is a small price to
60  // pay for more precise update below.
61  const Fractional leaving_squared_norm =
62  PreciseSquaredNorm(TransposedView(unit_row_left_inverse));
63  const Fractional old_squared_norm = edge_squared_norms_[leaving_row];
64  const Fractional estimated_edge_norms_accuracy =
65  (sqrt(leaving_squared_norm) - sqrt(old_squared_norm)) /
66  sqrt(leaving_squared_norm);
67  stats_.edge_norms_accuracy.Add(estimated_edge_norms_accuracy);
68  if (std::abs(estimated_edge_norms_accuracy) >
69  parameters_.recompute_edges_norm_threshold()) {
70  VLOG(1) << "Recomputing edge norms: " << sqrt(leaving_squared_norm)
71  << " vs " << sqrt(old_squared_norm);
72  recompute_edge_squared_norms_ = true;
73  return;
74  }
75 
76  const Fractional pivot = direction[leaving_row];
77  const Fractional new_leaving_squared_norm =
78  leaving_squared_norm / Square(pivot);
79 
80  // Update the norm.
81  int stat_lower_bounded_norms = 0;
82  for (const auto e : direction) {
83  // Note that the update formula used is important to maximize the precision.
84  // See Koberstein's PhD section 8.2.2.1.
85  edge_squared_norms_[e.row()] +=
86  e.coefficient() * (e.coefficient() * new_leaving_squared_norm -
87  2.0 / pivot * tau[e.row()]);
88 
89  // Avoid 0.0 norms (The 1e-4 is the value used by Koberstein).
90  // TODO(user): use a more precise lower bound depending on the column norm?
91  // We can do that with Cauchy-Swartz inequality:
92  // (edge . leaving_column)^2 = 1.0 < ||edge||^2 * ||leaving_column||^2
93  const Fractional kLowerBound = 1e-4;
94  if (edge_squared_norms_[e.row()] < kLowerBound) {
95  if (e.row() == leaving_row) continue;
96  edge_squared_norms_[e.row()] = kLowerBound;
97  ++stat_lower_bounded_norms;
98  }
99  }
100  edge_squared_norms_[leaving_row] = new_leaving_squared_norm;
101  IF_STATS_ENABLED(stats_.lower_bounded_norms.Add(stat_lower_bounded_norms));
102 }
103 
104 void DualEdgeNorms::ComputeEdgeSquaredNorms() {
105  SCOPED_TIME_STAT(&stats_);
106 
107  // Since we will do a lot of inversions, it is better to be as efficient and
108  // precise as possible by having a refactorized basis.
109  DCHECK(basis_factorization_.IsRefactorized());
110  const RowIndex num_rows = basis_factorization_.GetNumberOfRows();
111  edge_squared_norms_.resize(num_rows, 0.0);
112  for (RowIndex row(0); row < num_rows; ++row) {
113  edge_squared_norms_[row] = basis_factorization_.DualEdgeSquaredNorm(row);
114  }
115  recompute_edge_squared_norms_ = false;
116 }
117 
118 const DenseColumn& DualEdgeNorms::ComputeTau(
119  const ScatteredColumn& unit_row_left_inverse) {
120  SCOPED_TIME_STAT(&stats_);
121  const DenseColumn& result =
122  basis_factorization_.RightSolveForTau(unit_row_left_inverse);
123  IF_STATS_ENABLED(stats_.tau_density.Add(Density(Transpose(result))));
124  return result;
125 }
126 
127 } // namespace glop
128 } // namespace operations_research
#define DCHECK(condition)
Definition: base/logging.h:884
#define VLOG(verboselevel)
Definition: base/logging.h:978
const DenseColumn & RightSolveForTau(const ScatteredColumn &a) const
Fractional DualEdgeSquaredNorm(RowIndex row) const
void UpdateBeforeBasisPivot(ColIndex entering_col, RowIndex leaving_row, const ScatteredColumn &direction, const ScatteredRow &unit_row_left_inverse)
void UpdateDataOnBasisPermutation(const ColumnPermutation &col_perm)
DualEdgeNorms(const BasisFactorization &basis_factorization)
RowIndex row
Definition: markowitz.cc:175
Fractional Square(Fractional f)
Fractional PreciseSquaredNorm(const SparseColumn &v)
double Density(const DenseRow &row)
const DenseRow & Transpose(const DenseColumn &col)
void ApplyColumnPermutationToRowIndexedVector(const Permutation< ColIndex > &col_perm, RowIndexedVector *v)
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