20 #include "absl/strings/str_format.h"
32 : matrix_(nullptr), row_scale_(), col_scale_() {}
49 return row < row_scale_.
size() ? row_scale_[
row] : 1.0;
54 return col < col_scale_.
size() ? col_scale_[
col] : 1.0;
65 std::string SparseMatrixScaler::DebugInformationString()
const {
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"
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()));
97 DCHECK(matrix_ !=
nullptr);
101 if (min_magnitude == 0.0) {
105 VLOG(1) <<
"Before scaling:\n" << DebugInformationString();
106 if (method == GlopParameters::LINEAR_PROGRAM) {
109 if (lp_status.
ok()) {
117 const Fractional dynamic_range = max_magnitude / min_magnitude;
118 const Fractional kMaxDynamicRangeForGeometricScaling = 1e20;
119 if (dynamic_range < kMaxDynamicRangeForGeometricScaling) {
120 const int kScalingIterations = 4;
122 for (
int iteration = 0; iteration < kScalingIterations; ++iteration) {
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)) {
138 VLOG(1) <<
"Equilibration step: Rows scaled = " << rows_equilibrated
139 <<
", columns scaled = " << cols_equilibrated <<
"\n";
140 VLOG(1) << DebugInformationString();
150 for (I i(0); i < size; ++i) {
151 (*vector_to_scale)[i] *= scale[i];
154 for (I i(0); i < size; ++i) {
155 (*vector_to_scale)[i] /= scale[i];
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();
167 return (*scale_var_indices)[num];
172 DCHECK(row_vector !=
nullptr);
173 ScaleVector(col_scale_, up, row_vector);
178 DCHECK(column_vector !=
nullptr);
179 ScaleVector(row_scale_, up, column_vector);
183 DCHECK(matrix_ !=
nullptr);
187 const ColIndex num_cols = matrix_->
num_cols();
188 for (ColIndex
col(0);
col < num_cols; ++
col) {
190 const Fractional magnitude = fabs(e.coefficient());
191 if (magnitude != 0.0) {
192 sigma_square += magnitude * magnitude;
193 sigma_abs += magnitude;
198 if (n == 0.0)
return 0.0;
204 return (sigma_square - sigma_abs * sigma_abs / n) / n;
213 DCHECK(matrix_ !=
nullptr);
216 const ColIndex num_cols = matrix_->
num_cols();
217 for (ColIndex
col(0);
col < num_cols; ++
col) {
219 const Fractional magnitude = fabs(e.coefficient());
220 const RowIndex
row = e.row();
221 if (magnitude != 0.0) {
227 const RowIndex num_rows = matrix_->
num_rows();
229 for (RowIndex
row(0);
row < num_rows; ++
row) {
230 if (max_in_row[
row] == 0.0) {
231 scaling_factor[
row] = 1.0;
234 scaling_factor[
row] = sqrt(max_in_row[
row] * min_in_row[
row]);
237 return ScaleMatrixRows(scaling_factor);
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) {
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);
254 if (max_in_col != 0.0) {
256 ScaleMatrixColumn(
col, factor);
260 return num_cols_scaled;
269 DCHECK(matrix_ !=
nullptr);
270 const RowIndex num_rows = matrix_->
num_rows();
272 const ColIndex num_cols = matrix_->
num_cols();
273 for (ColIndex
col(0);
col < num_cols; ++
col) {
275 const Fractional magnitude = fabs(e.coefficient());
276 if (magnitude != 0.0) {
277 const RowIndex
row = e.row();
282 for (RowIndex
row(0);
row < num_rows; ++
row) {
283 if (max_magnitude[
row] == 0.0) {
284 max_magnitude[
row] = 1.0;
287 return ScaleMatrixRows(max_magnitude);
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) {
296 if (max_magnitude != 0.0) {
297 ScaleMatrixColumn(
col, max_magnitude);
301 return num_cols_scaled;
304 RowIndex SparseMatrixScaler::ScaleMatrixRows(
const DenseColumn& factors) {
306 DCHECK(matrix_ !=
nullptr);
307 const RowIndex num_rows = matrix_->
num_rows();
309 RowIndex num_rows_scaled(0);
310 for (RowIndex
row(0);
row < num_rows; ++
row) {
315 row_scale_[
row] *= factor;
319 const ColIndex num_cols = matrix_->
num_cols();
320 for (ColIndex
col(0);
col < num_cols; ++
col) {
322 if (column !=
nullptr) {
327 return num_rows_scaled;
330 void SparseMatrixScaler::ScaleMatrixColumn(ColIndex
col,
Fractional factor) {
332 DCHECK(matrix_ !=
nullptr);
333 col_scale_[
col] *= factor;
337 if (column !=
nullptr) {
344 DCHECK(matrix_ !=
nullptr);
345 const ColIndex num_cols = matrix_->
num_cols();
346 for (ColIndex
col(0);
col < num_cols; ++
col) {
351 if (column !=
nullptr) {
359 DCHECK(matrix_ !=
nullptr);
361 auto linear_program = absl::make_unique<LinearProgram>();
362 GlopParameters params;
363 auto simplex = absl::make_unique<RevisedSimplex>();
364 simplex->SetParameters(params);
388 const ColIndex beta = linear_program->CreateNewVariable();
391 linear_program->SetObjectiveCoefficient(beta, 1);
393 const ColIndex num_cols = matrix_->
num_cols();
394 for (ColIndex
col(0);
col < num_cols; ++
col) {
397 const ColIndex column_scale = CreateOrGetScaleIndex<ColIndex>(
398 col, linear_program.get(), &col_scale_var_indices);
405 const ColIndex row_scale = CreateOrGetScaleIndex<RowIndex>(
406 row, linear_program.get(), &row_scale_var_indices);
417 const RowIndex positive_constraint =
418 linear_program->CreateNewConstraint();
420 linear_program->SetConstraintBounds(positive_constraint, -log_magnitude,
423 linear_program->SetCoefficient(positive_constraint, row_scale, 1);
425 linear_program->SetCoefficient(positive_constraint, column_scale, 1);
427 linear_program->SetCoefficient(positive_constraint, beta, 1);
430 const RowIndex negative_constraint =
431 linear_program->CreateNewConstraint();
433 linear_program->SetConstraintBounds(negative_constraint, -
kInfinity,
436 linear_program->SetCoefficient(negative_constraint, row_scale, 1);
438 linear_program->SetCoefficient(negative_constraint, column_scale, 1);
440 linear_program->SetCoefficient(negative_constraint, beta, -1);
445 linear_program->AddSlackVariablesWhereNecessary(
false);
446 const Status simplex_status =
448 if (!simplex_status.
ok()) {
449 return simplex_status;
455 const ColIndex num_cols = matrix_->
num_cols();
456 for (ColIndex
col(0);
col < num_cols; ++
col) {
458 exp2(-simplex->GetVariableValue(CreateOrGetScaleIndex<ColIndex>(
459 col, linear_program.get(), &col_scale_var_indices)));
460 ScaleMatrixColumn(
col, column_scale);
462 const RowIndex num_rows = matrix_->
num_rows();
464 for (RowIndex
row(0);
row < num_rows; ++
row) {
466 exp2(-simplex->GetVariableValue(CreateOrGetScaleIndex<RowIndex>(
467 row, linear_program.get(), &row_scale_var_indices)));
469 ScaleMatrixRows(row_scale);
#define DCHECK_NE(val1, val2)
#define DCHECK_GE(val1, val2)
#define DCHECK(condition)
#define DCHECK_EQ(val1, val2)
#define VLOG(verboselevel)
void resize(size_type new_size)
static std::unique_ptr< TimeLimit > Infinite()
Creates a time limit object that uses infinite time for wall time, deterministic time and instruction...
Fractional EntryCoefficient(EntryIndex i) const
RowIndex EntryRow(EntryIndex i) const
SparseColumn * mutable_column(ColIndex col)
ColIndex num_cols() const
void ComputeMinAndMaxMagnitudes(Fractional *min_magnitude, Fractional *max_magnitude) const
RowIndex num_rows() const
const SparseColumn & column(ColIndex col) const
Fractional RowScalingFactor(RowIndex row) const
Fractional ColScalingFactor(ColIndex col) const
void ScaleColumnVector(bool up, DenseColumn *column_vector) const
void Init(SparseMatrix *matrix)
void ScaleRowVector(bool up, DenseRow *row_vector) const
ColIndex EquilibrateColumns()
RowIndex EquilibrateRows()
ColIndex ScaleColumnsGeometrically()
void Scale(GlopParameters::ScalingAlgorithm method)
Fractional ColUnscalingFactor(ColIndex col) const
Fractional RowUnscalingFactor(RowIndex row) const
RowIndex ScaleRowsGeometrically()
Fractional VarianceOfAbsoluteValueOfNonZeros() const
void ComponentWiseMultiply(const DenseVector &factors)
::util::IntegerRange< EntryIndex > AllEntryIndices() const
void DivideByConstant(Fractional factor)
void MultiplyByConstant(Fractional factor)
void ComponentWiseDivide(const DenseVector &factors)
typename Iterator::Entry Entry
const std::string & error_message() const
void resize(IntType size)
Fractional InfinityNorm(const DenseColumn &v)
Index ColToIntIndex(ColIndex col)
Index RowToIntIndex(RowIndex row)
const ColIndex kInvalidCol(-1)
static double ToDouble(double f)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
#define RETURN_IF_NULL(x)