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
18namespace operations_research {
19namespace glop {
20
22 : basis_factorization_(basis_factorization),
23 recompute_edge_squared_norms_(true) {}
24
26 return recompute_edge_squared_norms_;
27}
28
29void DualEdgeNorms::Clear() { recompute_edge_squared_norms_ = true; }
30
31void 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
104void 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
118const 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 PreciseSquaredNorm(const SparseColumn &v)
Fractional Square(Fractional f)
double Density(const DenseRow &row)
void ApplyColumnPermutationToRowIndexedVector(const Permutation< ColIndex > &col_perm, RowIndexedVector *v)
const DenseRow & Transpose(const DenseColumn &col)
const ScatteredRow & TransposedView(const ScatteredColumn &c)
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