OR-Tools  8.2
rank_one_update.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_RANK_ONE_UPDATE_H_
15#define OR_TOOLS_GLOP_RANK_ONE_UPDATE_H_
16
22
23namespace operations_research {
24namespace glop {
25
26// This class holds a matrix of the form T = I + u.Tr(v) where I is the
27// identity matrix and u and v are two column vectors of the same size as I. It
28// allows for efficient left and right solves with T. When T is non-singular,
29// it is easy to show that T^{-1} = I - 1 / mu * u.Tr(v) where
30// mu = 1.0 + Tr(v).u
31//
32// Note that when v is a unit vector, T is a regular Eta matrix and when u
33// is a unit vector, T is a row-wise Eta matrix.
34//
35// This is based on section 3.1 of:
36// Qi Huangfu, J. A. Julian Hall, "Novel update techniques for the revised
37// simplex method", 28 january 2013, Technical Report ERGO-13-0001
39 public:
40 // Rather than copying the vectors u and v, RankOneUpdateElementaryMatrix
41 // takes two columns of a provided CompactSparseMatrix which is used for
42 // storage. This has a couple of advantages, especially in the context of the
43 // RankOneUpdateFactorization below:
44 // - It uses less overall memory (and avoid allocation overhead).
45 // - It has a better cache behavior for the RankOneUpdateFactorization solves.
47 ColIndex u_index, ColIndex v_index,
48 Fractional u_dot_v)
49 : storage_(storage),
50 u_index_(u_index),
51 v_index_(v_index),
52 mu_(1.0 + u_dot_v) {}
53
54 // Returns whether or not this matrix is singular.
55 // Note that the RightSolve() and LeftSolve() function will fail if this is
56 // the case.
57 bool IsSingular() const { return mu_ == 0.0; }
58
59 // Solves T.x = rhs with rhs initialy in x (a column vector).
60 // The non-zeros version keeps track of the new non-zeros.
61 void RightSolve(DenseColumn* x) const {
63 const Fractional multiplier =
64 -storage_->ColumnScalarProduct(v_index_, Transpose(*x)) / mu_;
65 storage_->ColumnAddMultipleToDenseColumn(u_index_, multiplier, x);
66 }
69 const Fractional multiplier =
70 -storage_->ColumnScalarProduct(v_index_, Transpose(x->values)) / mu_;
71 if (multiplier != 0.0) {
72 storage_->ColumnAddMultipleToSparseScatteredColumn(u_index_, multiplier,
73 x);
74 }
75 }
76
77 // Solves y.T = rhs with rhs initialy in y (a row vector).
78 // The non-zeros version keeps track of the new non-zeros.
79 void LeftSolve(DenseRow* y) const {
81 const Fractional multiplier =
82 -storage_->ColumnScalarProduct(u_index_, *y) / mu_;
83 storage_->ColumnAddMultipleToDenseColumn(v_index_, multiplier,
84 reinterpret_cast<DenseColumn*>(y));
85 }
88 const Fractional multiplier =
89 -storage_->ColumnScalarProduct(u_index_, y->values) / mu_;
90 if (multiplier != 0.0) {
92 v_index_, multiplier, reinterpret_cast<ScatteredColumn*>(y));
93 }
94 }
95
96 // Computes T.x for a given column vector.
97 void RightMultiply(DenseColumn* x) const {
98 const Fractional multiplier =
99 storage_->ColumnScalarProduct(v_index_, Transpose(*x));
100 storage_->ColumnAddMultipleToDenseColumn(u_index_, multiplier, x);
101 }
102
103 // Computes y.T for a given row vector.
104 void LeftMultiply(DenseRow* y) const {
105 const Fractional multiplier = storage_->ColumnScalarProduct(u_index_, *y);
106 storage_->ColumnAddMultipleToDenseColumn(v_index_, multiplier,
107 reinterpret_cast<DenseColumn*>(y));
108 }
109
110 EntryIndex num_entries() const {
111 return storage_->column(u_index_).num_entries() +
112 storage_->column(v_index_).num_entries();
113 }
114
115 private:
116 // This is only used in debug mode.
117 Fractional ComputeUScalarV() const {
118 DenseColumn dense_u;
119 storage_->ColumnCopyToDenseColumn(u_index_, &dense_u);
120 return storage_->ColumnScalarProduct(v_index_, Transpose(dense_u));
121 }
122
123 // Note that we allow copy and assignment so we can store a
124 // RankOneUpdateElementaryMatrix in an STL container.
125 const CompactSparseMatrix* storage_;
126 ColIndex u_index_;
127 ColIndex v_index_;
128 Fractional mu_;
129};
130
131// A rank one update factorization corresponds to the product of k rank one
132// update elementary matrices, i.e. T = T_0.T_1. ... .T_{k-1}
134 public:
135 // TODO(user): make the 5% a parameter and share it between all the places
136 // that switch between a sparse/dense version.
137 RankOneUpdateFactorization() : hypersparse_ratio_(0.05) {}
138
139 // This is currently only visible for testing.
140 void set_hypersparse_ratio(double value) { hypersparse_ratio_ = value; }
141
142 // Deletes all elementary matrices of this factorization.
143 void Clear() {
144 elementary_matrices_.clear();
145 num_entries_ = 0;
146 }
147
148 // Updates the factorization.
149 void Update(const RankOneUpdateElementaryMatrix& update_matrix) {
150 elementary_matrices_.push_back(update_matrix);
151 num_entries_ += update_matrix.num_entries();
152 }
153
154 // Left-solves all systems from right to left, i.e. y_i = y_{i+1}.(T_i)^{-1}
155 void LeftSolve(DenseRow* y) const {
157 for (int i = elementary_matrices_.size() - 1; i >= 0; --i) {
158 elementary_matrices_[i].LeftSolve(y);
159 }
160 }
161
162 // Same as LeftSolve(), but if the given non_zeros are not empty, then all
163 // the new non-zeros in the result are appended to it.
166 if (y->non_zeros.empty()) {
167 LeftSolve(&y->values);
168 return;
169 }
170
171 // y->is_non_zero is always all false before and after this code.
174 bool use_dense = y->ShouldUseDenseIteration(hypersparse_ratio_);
175 for (int i = elementary_matrices_.size() - 1; i >= 0; --i) {
176 if (use_dense) {
177 elementary_matrices_[i].LeftSolve(&y->values);
178 } else {
179 elementary_matrices_[i].LeftSolveWithNonZeros(y);
180 use_dense = y->ShouldUseDenseIteration(hypersparse_ratio_);
181 }
182 }
183 y->ClearSparseMask();
184 y->ClearNonZerosIfTooDense(hypersparse_ratio_);
185 }
186
187 // Right-solves all systems from left to right, i.e. T_i.d_{i+1} = d_i
188 void RightSolve(DenseColumn* d) const {
190 const size_t end = elementary_matrices_.size();
191 for (int i = 0; i < end; ++i) {
192 elementary_matrices_[i].RightSolve(d);
193 }
194 }
195
196 // Same as RightSolve(), but if the given non_zeros are not empty, then all
197 // the new non-zeros in the result are appended to it.
200 if (d->non_zeros.empty()) {
201 RightSolve(&d->values);
202 return;
203 }
204
205 // d->is_non_zero is always all false before and after this code.
208 bool use_dense = d->ShouldUseDenseIteration(hypersparse_ratio_);
209 const size_t end = elementary_matrices_.size();
210 for (int i = 0; i < end; ++i) {
211 if (use_dense) {
212 elementary_matrices_[i].RightSolve(&d->values);
213 } else {
214 elementary_matrices_[i].RightSolveWithNonZeros(d);
215 use_dense = d->ShouldUseDenseIteration(hypersparse_ratio_);
216 }
217 }
218 d->ClearSparseMask();
219 d->ClearNonZerosIfTooDense(hypersparse_ratio_);
220 }
221
222 EntryIndex num_entries() const { return num_entries_; }
223
224 private:
225 double hypersparse_ratio_;
226 EntryIndex num_entries_;
227 std::vector<RankOneUpdateElementaryMatrix> elementary_matrices_;
228 DISALLOW_COPY_AND_ASSIGN(RankOneUpdateFactorization);
229};
230
231} // namespace glop
232} // namespace operations_research
233
234#endif // OR_TOOLS_GLOP_RANK_ONE_UPDATE_H_
#define DCHECK(condition)
Definition: base/logging.h:884
void ColumnCopyToDenseColumn(ColIndex col, DenseColumn *dense_column) const
Definition: sparse.h:418
void ColumnAddMultipleToSparseScatteredColumn(ColIndex col, Fractional multiplier, ScatteredColumn *column) const
Definition: sparse.h:405
Fractional ColumnScalarProduct(ColIndex col, const DenseRow &vector) const
Definition: sparse.h:382
void ColumnAddMultipleToDenseColumn(ColIndex col, Fractional multiplier, DenseColumn *dense_column) const
Definition: sparse.h:393
ColumnView column(ColIndex col) const
Definition: sparse.h:364
RankOneUpdateElementaryMatrix(const CompactSparseMatrix *storage, ColIndex u_index, ColIndex v_index, Fractional u_dot_v)
void Update(const RankOneUpdateElementaryMatrix &update_matrix)
int64 value
const DenseRow & Transpose(const DenseColumn &col)
bool IsAllFalse(const BoolVector &v)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
#define RETURN_IF_NULL(x)
Definition: return_macros.h:20
bool ShouldUseDenseIteration(double ratio_for_using_dense_representation) const
void ClearNonZerosIfTooDense(double ratio_for_using_dense_representation)
StrictITIVector< Index, bool > is_non_zero
StrictITIVector< Index, Fractional > values