23#include "absl/log/check.h"
24#include "absl/strings/str_format.h"
27#include "ortools/glop/parameters.pb.h"
42 : matrix_(nullptr), row_scales_(), col_scales_() {}
45 DCHECK(matrix !=
nullptr);
47 row_scales_.resize(matrix_->num_rows(), 1.0);
48 col_scales_.resize(matrix_->num_cols(), 1.0);
59 return row < row_scales_.size() ? row_scales_[row] : 1.0;
64 return col < col_scales_.size() ? col_scales_[col] : 1.0;
75std::string SparseMatrixScaler::DebugInformationString()
const {
79 DCHECK(!row_scales_.empty());
80 DCHECK(!col_scales_.empty());
84 const Fractional dynamic_range = max_magnitude / min_magnitude;
85 std::string output = absl::StrFormat(
86 "Min magnitude = %g, max magnitude = %g\n"
87 "Dynamic range = %g\n"
89 "Minimum row scale = %g, maximum row scale = %g\n"
90 "Minimum col scale = %g, maximum col scale = %g\n",
91 min_magnitude, max_magnitude, dynamic_range,
92 VarianceOfAbsoluteValueOfNonZeros(),
93 *std::min_element(row_scales_.begin(), row_scales_.end()),
94 *std::max_element(row_scales_.begin(), row_scales_.end()),
95 *std::min_element(col_scales_.begin(), col_scales_.end()),
96 *std::max_element(col_scales_.begin(), col_scales_.end()));
107 DCHECK(matrix_ !=
nullptr);
110 matrix_->ComputeMinAndMaxMagnitudes(&min_magnitude, &max_magnitude);
111 if (min_magnitude == 0.0) {
112 DCHECK_EQ(0.0, max_magnitude);
115 VLOG(1) <<
"Before scaling:\n" << DebugInformationString();
116 if (method == GlopParameters::LINEAR_PROGRAM) {
117 Status lp_status = LPScale();
119 if (lp_status.
ok()) {
122 VLOG(1) <<
"Error with LP scaling: " << lp_status.
error_message();
127 const Fractional dynamic_range = max_magnitude / min_magnitude;
128 const Fractional kMaxDynamicRangeForGeometricScaling = 1e20;
129 if (dynamic_range < kMaxDynamicRangeForGeometricScaling) {
130 const int kScalingIterations = 4;
132 for (
int iteration = 0; iteration < kScalingIterations; ++iteration) {
133 const RowIndex num_rows_scaled = ScaleRowsGeometrically();
134 const ColIndex num_cols_scaled = ScaleColumnsGeometrically();
135 const Fractional variance = VarianceOfAbsoluteValueOfNonZeros();
136 VLOG(1) <<
"Geometric scaling iteration " << iteration
137 <<
". Rows scaled = " << num_rows_scaled
138 <<
", columns scaled = " << num_cols_scaled <<
"\n";
139 VLOG(1) << DebugInformationString();
140 if (variance < kVarianceThreshold ||
141 (num_cols_scaled == 0 && num_rows_scaled == 0)) {
146 const RowIndex rows_equilibrated = EquilibrateRows();
147 const ColIndex cols_equilibrated = EquilibrateColumns();
148 VLOG(1) <<
"Equilibration step: Rows scaled = " << rows_equilibrated
149 <<
", columns scaled = " << cols_equilibrated <<
"\n";
150 VLOG(1) << DebugInformationString();
158 const I size(std::min(scale.size(), vector_to_scale->size()));
160 for (I i(0); i < size; ++i) {
161 (*vector_to_scale)[i] *= scale[i];
164 for (I i(0); i < size; ++i) {
165 (*vector_to_scale)[i] /= scale[i];
170template <
typename InputIndexType>
171ColIndex CreateOrGetScaleIndex(
173 util_intops::StrongVector<InputIndexType, ColIndex>* scale_var_indices) {
174 if ((*scale_var_indices)[num] == -1) {
175 (*scale_var_indices)[num] = lp->CreateNewVariable();
177 return (*scale_var_indices)[num];
182 DCHECK(row_vector !=
nullptr);
183 ScaleVector(col_scales_, up, row_vector);
188 DCHECK(column_vector !=
nullptr);
189 ScaleVector(row_scales_, up, column_vector);
192Fractional SparseMatrixScaler::VarianceOfAbsoluteValueOfNonZeros()
const {
193 DCHECK(matrix_ !=
nullptr);
197 const ColIndex num_cols = matrix_->
num_cols();
198 for (ColIndex col(0); col < num_cols; ++col) {
199 for (
const SparseColumn::Entry e : matrix_->
column(col)) {
200 const Fractional magnitude = fabs(e.coefficient());
201 sigma_square += magnitude * magnitude;
202 sigma_abs += magnitude;
206 if (n == 0.0)
return 0.0;
213 return (sigma_square - sigma_abs * sigma_abs / n) / n;
221RowIndex SparseMatrixScaler::ScaleRowsGeometrically() {
222 DCHECK(matrix_ !=
nullptr);
225 const ColIndex num_cols = matrix_->num_cols();
226 for (ColIndex col(0); col < num_cols; ++col) {
227 for (
const SparseColumn::Entry e : matrix_->column(col)) {
228 const Fractional magnitude = fabs(e.coefficient());
229 const RowIndex row = e.row();
230 if (magnitude != 0.0) {
231 max_in_row[row] = std::max(max_in_row[row], magnitude);
232 min_in_row[row] = std::min(min_in_row[row], magnitude);
236 const RowIndex num_rows = matrix_->num_rows();
238 for (RowIndex row(0); row < num_rows; ++row) {
239 if (max_in_row[row] == 0.0) {
240 scaling_factor[row] = 1.0;
243 scaling_factor[row] = std::sqrt(max_in_row[row] * min_in_row[row]);
246 return ScaleMatrixRows(scaling_factor);
249ColIndex SparseMatrixScaler::ScaleColumnsGeometrically() {
250 DCHECK(matrix_ !=
nullptr);
251 ColIndex num_cols_scaled(0);
252 const ColIndex num_cols = matrix_->num_cols();
253 for (ColIndex col(0); col < num_cols; ++col) {
256 for (
const SparseColumn::Entry e : matrix_->column(col)) {
257 const Fractional magnitude = fabs(e.coefficient());
258 if (magnitude != 0.0) {
259 max_in_col = std::max(max_in_col, magnitude);
260 min_in_col = std::min(min_in_col, magnitude);
263 if (max_in_col != 0.0) {
265 ScaleMatrixColumn(col, factor);
269 return num_cols_scaled;
277RowIndex SparseMatrixScaler::EquilibrateRows() {
278 DCHECK(matrix_ !=
nullptr);
279 const RowIndex num_rows = matrix_->num_rows();
281 const ColIndex num_cols = matrix_->num_cols();
282 for (ColIndex col(0); col < num_cols; ++col) {
283 for (
const SparseColumn::Entry e : matrix_->column(col)) {
284 const Fractional magnitude = fabs(e.coefficient());
285 if (magnitude != 0.0) {
286 const RowIndex row = e.row();
287 max_magnitudes[row] = std::max(max_magnitudes[row], magnitude);
291 for (RowIndex row(0); row < num_rows; ++row) {
292 if (max_magnitudes[row] == 0.0) {
293 max_magnitudes[row] = 1.0;
296 return ScaleMatrixRows(max_magnitudes);
299ColIndex SparseMatrixScaler::EquilibrateColumns() {
300 DCHECK(matrix_ !=
nullptr);
301 ColIndex num_cols_scaled(0);
302 const ColIndex num_cols = matrix_->num_cols();
303 for (ColIndex col(0); col < num_cols; ++col) {
305 if (max_magnitude != 0.0 && max_magnitude != 1.0) {
306 ScaleMatrixColumn(col, max_magnitude);
310 return num_cols_scaled;
313RowIndex SparseMatrixScaler::ScaleMatrixRows(
const DenseColumn& factors) {
315 DCHECK(matrix_ !=
nullptr);
316 const RowIndex num_rows = matrix_->num_rows();
317 DCHECK_EQ(num_rows, factors.size());
318 RowIndex num_rows_scaled(0);
319 for (RowIndex row(0); row < num_rows; ++row) {
321 DCHECK_NE(0.0, factor);
324 row_scales_[row] *= factor;
328 const ColIndex num_cols = matrix_->num_cols();
329 for (ColIndex col(0); col < num_cols; ++col) {
330 matrix_->mutable_column(col)->ComponentWiseDivide(factors);
333 return num_rows_scaled;
336void SparseMatrixScaler::ScaleMatrixColumn(ColIndex col,
Fractional factor) {
338 DCHECK(matrix_ !=
nullptr);
339 DCHECK_NE(0.0, factor);
340 col_scales_[col] *= factor;
341 matrix_->mutable_column(col)->DivideByConstant(factor);
344Status SparseMatrixScaler::LPScale() {
345 DCHECK(matrix_ !=
nullptr);
347 auto linear_program = std::make_unique<LinearProgram>();
348 GlopParameters params;
349 auto simplex = std::make_unique<RevisedSimplex>();
350 simplex->SetParameters(params);
370 util_intops::StrongVector<ColIndex, ColIndex> col_scale_var_indices;
371 util_intops::StrongVector<RowIndex, ColIndex> row_scale_var_indices;
374 const ColIndex beta = linear_program->CreateNewVariable();
377 linear_program->SetObjectiveCoefficient(beta, 1);
379 const ColIndex num_cols = matrix_->num_cols();
380 for (ColIndex col(0); col < num_cols; ++col) {
382 const ColIndex column_scale = CreateOrGetScaleIndex<ColIndex>(
383 col, linear_program.get(), &col_scale_var_indices);
385 for (
const SparseColumn::Entry e : matrix_->column(col)) {
386 const Fractional log_magnitude = log2(std::abs(e.coefficient()));
387 const RowIndex row = e.row();
390 const ColIndex row_scale = CreateOrGetScaleIndex<RowIndex>(
391 row, linear_program.get(), &row_scale_var_indices);
404 const RowIndex positive_constraint =
405 linear_program->CreateNewConstraint();
407 linear_program->SetConstraintBounds(positive_constraint, -log_magnitude,
410 linear_program->SetCoefficient(positive_constraint, row_scale, 1);
412 linear_program->SetCoefficient(positive_constraint, column_scale, 1);
414 linear_program->SetCoefficient(positive_constraint, beta, 1);
417 const RowIndex negative_constraint =
418 linear_program->CreateNewConstraint();
420 linear_program->SetConstraintBounds(negative_constraint, -
kInfinity,
423 linear_program->SetCoefficient(negative_constraint, row_scale, 1);
425 linear_program->SetCoefficient(negative_constraint, column_scale, 1);
427 linear_program->SetCoefficient(negative_constraint, beta, -1);
432 linear_program->AddSlackVariablesWhereNecessary(
false);
433 const Status simplex_status =
435 if (!simplex_status.ok()) {
436 return simplex_status;
442 const ColIndex num_cols = matrix_->num_cols();
443 for (ColIndex col(0); col < num_cols; ++col) {
445 exp2(-simplex->GetVariableValue(CreateOrGetScaleIndex<ColIndex>(
446 col, linear_program.get(), &col_scale_var_indices)));
447 ScaleMatrixColumn(col, column_scale);
449 const RowIndex num_rows = matrix_->num_rows();
451 for (RowIndex row(0); row < num_rows; ++row) {
453 exp2(-simplex->GetVariableValue(CreateOrGetScaleIndex<RowIndex>(
454 row, linear_program.get(), &row_scale_var_indices)));
456 ScaleMatrixRows(row_scale);
static std::unique_ptr< TimeLimit > Infinite()
void ScaleColumnVector(bool up, DenseColumn *column_vector) const
void Init(SparseMatrix *matrix)
Fractional RowScalingFactor(RowIndex row) const
void ScaleRowVector(bool up, DenseRow *row_vector) const
Fractional ColScalingFactor(ColIndex col) const
Fractional RowUnscalingFactor(RowIndex row) const
void Scale(GlopParameters::ScalingAlgorithm method)
Scales the matrix.
Fractional ColUnscalingFactor(ColIndex col) const
void ComputeMinAndMaxMagnitudes(Fractional *min_magnitude, Fractional *max_magnitude) const
const SparseColumn & column(ColIndex col) const
Access the underlying sparse columns.
ColIndex num_cols() const
static Status OK()
Improves readability but identical to 0-arg constructor.
const std::string & error_message() const
void resize(size_type new_size)
constexpr ColIndex kInvalidCol(-1)
constexpr double kInfinity
Infinity for type Fractional.
Index ColToIntIndex(ColIndex col)
Get the integer index corresponding to the col.
Fractional InfinityNorm(const DenseColumn &v)
Returns the maximum of the coefficients of 'v'.
StrictITIVector< RowIndex, Fractional > DenseColumn
Column-vector types. Column-vector types are indexed by a row index.
StrictITIVector< ColIndex, Fractional > DenseRow
Row-vector types. Row-vector types are indexed by a column index.
Index RowToIntIndex(RowIndex row)
Get the integer index corresponding to the row.
static double ToDouble(double f)
In SWIG mode, we don't want anything besides these top-level includes.
#define RETURN_IF_NULL(x)