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"
27
28namespace operations_research {
29namespace 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
65std::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
90void 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
143namespace {
144template <class I>
145void 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
160template <typename InputIndexType>
161ColIndex 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
171void 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
304RowIndex 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
330void 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)
void DivideByConstant(Fractional factor)
void MultiplyByConstant(Fractional factor)
void ComponentWiseDivide(const DenseVector &factors)
::util::IntegerRange< EntryIndex > AllEntryIndices() const
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)
const ColIndex kInvalidCol(-1)
Index RowToIntIndex(RowIndex row)
Definition: lp_types.h:57
static double ToDouble(double f)
Definition: lp_types.h:68
Index ColToIntIndex(ColIndex col)
Definition: lp_types.h:54
const double kInfinity
Definition: lp_types.h:83
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