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 
17 #include "ortools/base/logging.h"
21 #include "ortools/lp_data/sparse.h"
22 
23 namespace operations_research {
24 namespace 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 {
62  DCHECK(!IsSingular());
63  const Fractional multiplier =
64  -storage_->ColumnScalarProduct(v_index_, Transpose(*x)) / mu_;
65  storage_->ColumnAddMultipleToDenseColumn(u_index_, multiplier, x);
66  }
68  DCHECK(!IsSingular());
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 {
80  DCHECK(!IsSingular());
81  const Fractional multiplier =
82  -storage_->ColumnScalarProduct(u_index_, *y) / mu_;
83  storage_->ColumnAddMultipleToDenseColumn(v_index_, multiplier,
84  reinterpret_cast<DenseColumn*>(y));
85  }
87  DCHECK(!IsSingular());
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 {
156  RETURN_IF_NULL(y);
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.
165  RETURN_IF_NULL(y);
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 {
189  RETURN_IF_NULL(d);
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.
199  RETURN_IF_NULL(d);
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
bool IsAllFalse(const BoolVector &v)
const DenseRow & Transpose(const DenseColumn &col)
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