OR-Tools  8.2
matrix_scaler.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 <algorithm>
17 #include <cmath>
18 #include <vector>
19 
20 #include "absl/strings/str_format.h"
21 #include "ortools/base/logging.h"
25 #include "ortools/lp_data/sparse.h"
27 
28 namespace operations_research {
29 namespace glop {
30 
32  : matrix_(nullptr), row_scale_(), col_scale_() {}
33 
35  DCHECK(matrix != nullptr);
36  matrix_ = matrix;
37  row_scale_.resize(matrix_->num_rows(), 1.0);
38  col_scale_.resize(matrix_->num_cols(), 1.0);
39 }
40 
42  matrix_ = nullptr;
43  row_scale_.clear();
44  col_scale_.clear();
45 }
46 
48  DCHECK_GE(row, 0);
49  return row < row_scale_.size() ? row_scale_[row] : 1.0;
50 }
51 
53  DCHECK_GE(col, 0);
54  return col < col_scale_.size() ? col_scale_[col] : 1.0;
55 }
56 
58  return 1.0 / RowUnscalingFactor(row);
59 }
60 
62  return 1.0 / ColUnscalingFactor(col);
63 }
64 
65 std::string SparseMatrixScaler::DebugInformationString() const {
66  // Note that some computations are redundant with the computations made in
67  // some callees, but we do not care as this function is supposed to be called
68  // with FLAGS_v set to 1.
69  DCHECK(!row_scale_.empty());
70  DCHECK(!col_scale_.empty());
71  Fractional max_magnitude;
72  Fractional min_magnitude;
73  matrix_->ComputeMinAndMaxMagnitudes(&min_magnitude, &max_magnitude);
74  const Fractional dynamic_range = max_magnitude / min_magnitude;
75  std::string output = absl::StrFormat(
76  "Min magnitude = %g, max magnitude = %g\n"
77  "Dynamic range = %g\n"
78  "Variance = %g\n"
79  "Minimum row scale = %g, maximum row scale = %g\n"
80  "Minimum col scale = %g, maximum col scale = %g\n",
81  min_magnitude, max_magnitude, dynamic_range,
83  *std::min_element(row_scale_.begin(), row_scale_.end()),
84  *std::max_element(row_scale_.begin(), row_scale_.end()),
85  *std::min_element(col_scale_.begin(), col_scale_.end()),
86  *std::max_element(col_scale_.begin(), col_scale_.end()));
87  return output;
88 }
89 
90 void SparseMatrixScaler::Scale(GlopParameters::ScalingAlgorithm method) {
91  // This is an implementation of the algorithm described in
92  // Benichou, M., Gauthier, J-M., Hentges, G., and Ribiere, G.,
93  // "The efficient solution of large-scale linear programming problems —
94  // some algorithmic techniques and computational results,"
95  // Mathematical Programming 13(3) (December 1977).
96  // http://www.springerlink.com/content/j3367676856m0064/
97  DCHECK(matrix_ != nullptr);
98  Fractional max_magnitude;
99  Fractional min_magnitude;
100  matrix_->ComputeMinAndMaxMagnitudes(&min_magnitude, &max_magnitude);
101  if (min_magnitude == 0.0) {
102  DCHECK_EQ(0.0, max_magnitude);
103  return; // Null matrix: nothing to do.
104  }
105  VLOG(1) << "Before scaling:\n" << DebugInformationString();
106  if (method == GlopParameters::LINEAR_PROGRAM) {
107  Status lp_status = LPScale();
108  // Revert to the default scaling method if there is an error with the LP.
109  if (lp_status.ok()) {
110  return;
111  } else {
112  VLOG(1) << "Error with LP scaling: " << lp_status.error_message();
113  }
114  }
115  // TODO(user): Decide precisely for which value of dynamic range we should cut
116  // off geometric scaling.
117  const Fractional dynamic_range = max_magnitude / min_magnitude;
118  const Fractional kMaxDynamicRangeForGeometricScaling = 1e20;
119  if (dynamic_range < kMaxDynamicRangeForGeometricScaling) {
120  const int kScalingIterations = 4;
121  const Fractional kVarianceThreshold(10.0);
122  for (int iteration = 0; iteration < kScalingIterations; ++iteration) {
123  const RowIndex num_rows_scaled = ScaleRowsGeometrically();
124  const ColIndex num_cols_scaled = ScaleColumnsGeometrically();
126  VLOG(1) << "Geometric scaling iteration " << iteration
127  << ". Rows scaled = " << num_rows_scaled
128  << ", columns scaled = " << num_cols_scaled << "\n";
129  VLOG(1) << DebugInformationString();
130  if (variance < kVarianceThreshold ||
131  (num_cols_scaled == 0 && num_rows_scaled == 0)) {
132  break;
133  }
134  }
135  }
136  RowIndex rows_equilibrated = EquilibrateRows();
137  ColIndex cols_equilibrated = EquilibrateColumns();
138  VLOG(1) << "Equilibration step: Rows scaled = " << rows_equilibrated
139  << ", columns scaled = " << cols_equilibrated << "\n";
140  VLOG(1) << DebugInformationString();
141 }
142 
143 namespace {
144 template <class I>
145 void ScaleVector(const absl::StrongVector<I, Fractional>& scale, bool up,
146  absl::StrongVector<I, Fractional>* vector_to_scale) {
147  RETURN_IF_NULL(vector_to_scale);
148  const I size(std::min(scale.size(), vector_to_scale->size()));
149  if (up) {
150  for (I i(0); i < size; ++i) {
151  (*vector_to_scale)[i] *= scale[i];
152  }
153  } else {
154  for (I i(0); i < size; ++i) {
155  (*vector_to_scale)[i] /= scale[i];
156  }
157  }
158 }
159 
160 template <typename InputIndexType>
161 ColIndex CreateOrGetScaleIndex(
162  InputIndexType num, LinearProgram* lp,
164  if ((*scale_var_indices)[num] == -1) {
165  (*scale_var_indices)[num] = lp->CreateNewVariable();
166  }
167  return (*scale_var_indices)[num];
168 }
169 } // anonymous namespace
170 
171 void SparseMatrixScaler::ScaleRowVector(bool up, DenseRow* row_vector) const {
172  DCHECK(row_vector != nullptr);
173  ScaleVector(col_scale_, up, row_vector);
174 }
175 
177  DenseColumn* column_vector) const {
178  DCHECK(column_vector != nullptr);
179  ScaleVector(row_scale_, up, column_vector);
180 }
181 
183  DCHECK(matrix_ != nullptr);
184  Fractional sigma_square(0.0);
185  Fractional sigma_abs(0.0);
186  double n = 0.0; // n is used in a calculation involving doubles.
187  const ColIndex num_cols = matrix_->num_cols();
188  for (ColIndex col(0); col < num_cols; ++col) {
189  for (const SparseColumn::Entry e : matrix_->column(col)) {
190  const Fractional magnitude = fabs(e.coefficient());
191  if (magnitude != 0.0) {
192  sigma_square += magnitude * magnitude;
193  sigma_abs += magnitude;
194  ++n;
195  }
196  }
197  }
198  if (n == 0.0) return 0.0;
199  // Since we know all the population (the non-zeros) and we are not using a
200  // sample, the variance is defined as below.
201  // For an explanation, see:
202  // http://en.wikipedia.org/wiki/Variance
203  // #Population_variance_and_sample_variance
204  return (sigma_square - sigma_abs * sigma_abs / n) / n;
205 }
206 
207 // For geometric scaling, we compute the maximum and minimum magnitudes
208 // of non-zeros in a row (resp. column). Let us denote these numbers as
209 // max and min. We then scale the row (resp. column) by dividing the
210 // coefficients by sqrt(min * max).
211 
213  DCHECK(matrix_ != nullptr);
214  DenseColumn max_in_row(matrix_->num_rows(), 0.0);
215  DenseColumn min_in_row(matrix_->num_rows(), kInfinity);
216  const ColIndex num_cols = matrix_->num_cols();
217  for (ColIndex col(0); col < num_cols; ++col) {
218  for (const SparseColumn::Entry e : matrix_->column(col)) {
219  const Fractional magnitude = fabs(e.coefficient());
220  const RowIndex row = e.row();
221  if (magnitude != 0.0) {
222  max_in_row[row] = std::max(max_in_row[row], magnitude);
223  min_in_row[row] = std::min(min_in_row[row], magnitude);
224  }
225  }
226  }
227  const RowIndex num_rows = matrix_->num_rows();
228  DenseColumn scaling_factor(num_rows, 0.0);
229  for (RowIndex row(0); row < num_rows; ++row) {
230  if (max_in_row[row] == 0.0) {
231  scaling_factor[row] = 1.0;
232  } else {
233  DCHECK_NE(kInfinity, min_in_row[row]);
234  scaling_factor[row] = sqrt(max_in_row[row] * min_in_row[row]);
235  }
236  }
237  return ScaleMatrixRows(scaling_factor);
238 }
239 
241  DCHECK(matrix_ != nullptr);
242  ColIndex num_cols_scaled(0);
243  const ColIndex num_cols = matrix_->num_cols();
244  for (ColIndex col(0); col < num_cols; ++col) {
245  Fractional max_in_col(0.0);
246  Fractional min_in_col(kInfinity);
247  for (const SparseColumn::Entry e : matrix_->column(col)) {
248  const Fractional magnitude = fabs(e.coefficient());
249  if (magnitude != 0.0) {
250  max_in_col = std::max(max_in_col, magnitude);
251  min_in_col = std::min(min_in_col, magnitude);
252  }
253  }
254  if (max_in_col != 0.0) {
255  const Fractional factor(sqrt(ToDouble(max_in_col * min_in_col)));
256  ScaleMatrixColumn(col, factor);
257  num_cols_scaled++;
258  }
259  }
260  return num_cols_scaled;
261 }
262 
263 // For equilibration, we compute the maximum magnitude of non-zeros
264 // in a row (resp. column), and then scale the row (resp. column) by dividing
265 // the coefficients this maximum magnitude.
266 // This brings the largest coefficient in a row equal to 1.0.
267 
269  DCHECK(matrix_ != nullptr);
270  const RowIndex num_rows = matrix_->num_rows();
271  DenseColumn max_magnitude(num_rows, 0.0);
272  const ColIndex num_cols = matrix_->num_cols();
273  for (ColIndex col(0); col < num_cols; ++col) {
274  for (const SparseColumn::Entry e : matrix_->column(col)) {
275  const Fractional magnitude = fabs(e.coefficient());
276  if (magnitude != 0.0) {
277  const RowIndex row = e.row();
278  max_magnitude[row] = std::max(max_magnitude[row], magnitude);
279  }
280  }
281  }
282  for (RowIndex row(0); row < num_rows; ++row) {
283  if (max_magnitude[row] == 0.0) {
284  max_magnitude[row] = 1.0;
285  }
286  }
287  return ScaleMatrixRows(max_magnitude);
288 }
289 
291  DCHECK(matrix_ != nullptr);
292  ColIndex num_cols_scaled(0);
293  const ColIndex num_cols = matrix_->num_cols();
294  for (ColIndex col(0); col < num_cols; ++col) {
295  const Fractional max_magnitude = InfinityNorm(matrix_->column(col));
296  if (max_magnitude != 0.0) {
297  ScaleMatrixColumn(col, max_magnitude);
298  num_cols_scaled++;
299  }
300  }
301  return num_cols_scaled;
302 }
303 
304 RowIndex SparseMatrixScaler::ScaleMatrixRows(const DenseColumn& factors) {
305  // Matrix rows are scaled by dividing their coefficients by factors[row].
306  DCHECK(matrix_ != nullptr);
307  const RowIndex num_rows = matrix_->num_rows();
308  DCHECK_EQ(num_rows, factors.size());
309  RowIndex num_rows_scaled(0);
310  for (RowIndex row(0); row < num_rows; ++row) {
311  const Fractional factor = factors[row];
312  DCHECK_NE(0.0, factor);
313  if (factor != 1.0) {
314  ++num_rows_scaled;
315  row_scale_[row] *= factor;
316  }
317  }
318 
319  const ColIndex num_cols = matrix_->num_cols();
320  for (ColIndex col(0); col < num_cols; ++col) {
321  SparseColumn* const column = matrix_->mutable_column(col);
322  if (column != nullptr) {
323  column->ComponentWiseDivide(factors);
324  }
325  }
326 
327  return num_rows_scaled;
328 }
329 
330 void SparseMatrixScaler::ScaleMatrixColumn(ColIndex col, Fractional factor) {
331  // A column is scaled by dividing by factor.
332  DCHECK(matrix_ != nullptr);
333  col_scale_[col] *= factor;
334  DCHECK_NE(0.0, factor);
335 
336  SparseColumn* const column = matrix_->mutable_column(col);
337  if (column != nullptr) {
338  column->DivideByConstant(factor);
339  }
340 }
341 
343  // Unscaling is easier than scaling since all scaling factors are stored.
344  DCHECK(matrix_ != nullptr);
345  const ColIndex num_cols = matrix_->num_cols();
346  for (ColIndex col(0); col < num_cols; ++col) {
347  const Fractional column_scale = col_scale_[col];
348  DCHECK_NE(0.0, column_scale);
349 
350  SparseColumn* const column = matrix_->mutable_column(col);
351  if (column != nullptr) {
352  column->MultiplyByConstant(column_scale);
353  column->ComponentWiseMultiply(row_scale_);
354  }
355  }
356 }
357 
359  DCHECK(matrix_ != nullptr);
360 
361  auto linear_program = absl::make_unique<LinearProgram>();
362  GlopParameters params;
363  auto simplex = absl::make_unique<RevisedSimplex>();
364  simplex->SetParameters(params);
365 
366  // Begin linear program construction.
367  // Beta represents the largest distance from zero among the constraint pairs.
368  // It resembles a slack variable because the 'objective' of each constraint is
369  // to cancel out the log "w" of the original nonzero |a_ij| (a.k.a. |a_rc|).
370  // Approaching 0 by addition in log space is the same as approaching 1 by
371  // multiplication in linear space. Hence, each variable's log magnitude is
372  // subtracted from the log row scale and log column scale. If the sum is
373  // positive, the positive constraint is trivially satisfied, but the negative
374  // constraint will determine the minimum necessary value of beta for that
375  // variable and scaling factors, and vice versa.
376  // For an MxN matrix, the resulting scaling LP has M+N+1 variables and
377  // O(M*N) constraints (2*M*N at maximum). As a result, using this LP to scale
378  // another linear program, will typically increase the time to
379  // optimization by a factor of 4, and has increased the time of some benchmark
380  // LPs by up to 10.
381 
382  // Indices to variables in the LinearProgram populated by
383  // GenerateLinearProgram.
384  absl::StrongVector<ColIndex, ColIndex> col_scale_var_indices;
385  absl::StrongVector<RowIndex, ColIndex> row_scale_var_indices;
386  row_scale_var_indices.resize(RowToIntIndex(matrix_->num_rows()), kInvalidCol);
387  col_scale_var_indices.resize(ColToIntIndex(matrix_->num_cols()), kInvalidCol);
388  const ColIndex beta = linear_program->CreateNewVariable();
389  linear_program->SetVariableBounds(beta, -kInfinity, kInfinity);
390  // Default objective is to minimize.
391  linear_program->SetObjectiveCoefficient(beta, 1);
392  matrix_->CleanUp();
393  const ColIndex num_cols = matrix_->num_cols();
394  for (ColIndex col(0); col < num_cols; ++col) {
395  SparseColumn* const column = matrix_->mutable_column(col);
396  // This is the variable representing the log of the scale factor for col.
397  const ColIndex column_scale = CreateOrGetScaleIndex<ColIndex>(
398  col, linear_program.get(), &col_scale_var_indices);
399  linear_program->SetVariableBounds(column_scale, -kInfinity, kInfinity);
400  for (EntryIndex i : column->AllEntryIndices()) {
401  const Fractional log_magnitude =
402  log2(std::abs(column->EntryCoefficient(i)));
403  const RowIndex row = column->EntryRow(i);
404  // This is the variable representing the log of the scale factor for row.
405  const ColIndex row_scale = CreateOrGetScaleIndex<RowIndex>(
406  row, linear_program.get(), &row_scale_var_indices);
407 
408  linear_program->SetVariableBounds(row_scale, -kInfinity, kInfinity);
409  // This is derived from the formulation in
410  // min β
411  // Subject to:
412  // ∀ c∈C, v∈V, p_{c,v} ≠ 0.0, w_{c,v} + s^{var}_v + s^{comb}_c + β ≥ 0.0
413  // ∀ c∈C, v∈V, p_{c,v} ≠ 0.0, w_{c,v} + s^{var}_v + s^{comb}_c ≤ β
414  // If a variable is integer, its scale factor is zero.
415 
416  // Start with the constraint w_cv + s_c + s_v + beta >= 0.
417  const RowIndex positive_constraint =
418  linear_program->CreateNewConstraint();
419  // Subtract the constant w_cv from both sides.
420  linear_program->SetConstraintBounds(positive_constraint, -log_magnitude,
421  kInfinity);
422  // +s_c, meaning (log) scale of the constraint C, pointed by row_scale.
423  linear_program->SetCoefficient(positive_constraint, row_scale, 1);
424  // +s_v, meaning (log) scale of the variable V, pointed by column_scale.
425  linear_program->SetCoefficient(positive_constraint, column_scale, 1);
426  // +beta
427  linear_program->SetCoefficient(positive_constraint, beta, 1);
428 
429  // Construct the constraint w_cv + s_c + s_v <= beta.
430  const RowIndex negative_constraint =
431  linear_program->CreateNewConstraint();
432  // Subtract w (and beta) from both sides.
433  linear_program->SetConstraintBounds(negative_constraint, -kInfinity,
434  -log_magnitude);
435  // +s_c, meaning (log) scale of the constraint C, pointed by row_scale.
436  linear_program->SetCoefficient(negative_constraint, row_scale, 1);
437  // +s_v, meaning (log) scale of the variable V, pointed by column_scale.
438  linear_program->SetCoefficient(negative_constraint, column_scale, 1);
439  // -beta
440  linear_program->SetCoefficient(negative_constraint, beta, -1);
441  }
442  }
443  // End linear program construction.
444 
445  linear_program->AddSlackVariablesWhereNecessary(false);
446  const Status simplex_status =
447  simplex->Solve(*linear_program, TimeLimit::Infinite().get());
448  if (!simplex_status.ok()) {
449  return simplex_status;
450  } else {
451  // Now the solution variables can be interpreted and translated from log
452  // space.
453  // For each row scale, unlog it and scale the constraints and constraint
454  // bounds.
455  const ColIndex num_cols = matrix_->num_cols();
456  for (ColIndex col(0); col < num_cols; ++col) {
457  const Fractional column_scale =
458  exp2(-simplex->GetVariableValue(CreateOrGetScaleIndex<ColIndex>(
459  col, linear_program.get(), &col_scale_var_indices)));
460  ScaleMatrixColumn(col, column_scale);
461  }
462  const RowIndex num_rows = matrix_->num_rows();
463  DenseColumn row_scale(num_rows, 0.0);
464  for (RowIndex row(0); row < num_rows; ++row) {
465  row_scale[row] =
466  exp2(-simplex->GetVariableValue(CreateOrGetScaleIndex<RowIndex>(
467  row, linear_program.get(), &row_scale_var_indices)));
468  }
469  ScaleMatrixRows(row_scale);
470  return Status::OK();
471  }
472 }
473 
474 } // namespace glop
475 } // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define DCHECK(condition)
Definition: base/logging.h:884
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
#define VLOG(verboselevel)
Definition: base/logging.h:978
void resize(size_type new_size)
size_type size() const
bool empty() const
static std::unique_ptr< TimeLimit > Infinite()
Creates a time limit object that uses infinite time for wall time, deterministic time and instruction...
Definition: time_limit.h:134
Fractional EntryCoefficient(EntryIndex i) const
Definition: sparse_column.h:52
RowIndex EntryRow(EntryIndex i) const
Definition: sparse_column.h:51
SparseColumn * mutable_column(ColIndex col)
Definition: sparse.h:181
void ComputeMinAndMaxMagnitudes(Fractional *min_magnitude, Fractional *max_magnitude) const
Definition: sparse.cc:369
const SparseColumn & column(ColIndex col) const
Definition: sparse.h:180
Fractional RowScalingFactor(RowIndex row) const
Fractional ColScalingFactor(ColIndex col) const
void ScaleColumnVector(bool up, DenseColumn *column_vector) const
void ScaleRowVector(bool up, DenseRow *row_vector) const
void Scale(GlopParameters::ScalingAlgorithm method)
Fractional ColUnscalingFactor(ColIndex col) const
Fractional RowUnscalingFactor(RowIndex row) const
void ComponentWiseMultiply(const DenseVector &factors)
::util::IntegerRange< EntryIndex > AllEntryIndices() const
void DivideByConstant(Fractional factor)
void MultiplyByConstant(Fractional factor)
void ComponentWiseDivide(const DenseVector &factors)
static const Status OK()
Definition: status.h:54
const std::string & error_message() const
Definition: status.h:58
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
Fractional InfinityNorm(const DenseColumn &v)
Index ColToIntIndex(ColIndex col)
Definition: lp_types.h:54
Index RowToIntIndex(RowIndex row)
Definition: lp_types.h:57
const double kInfinity
Definition: lp_types.h:83
const ColIndex kInvalidCol(-1)
static double ToDouble(double f)
Definition: lp_types.h:68
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