OR-Tools  9.6
sharded_optimization_utils_test.cc
Go to the documentation of this file.
1 // Copyright 2010-2022 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 #include "ortools/pdlp/sharded_optimization_utils.h"
15 
16 #include <cmath>
17 #include <cstdint>
18 #include <optional>
19 #include <random>
20 #include <utility>
21 #include <vector>
22 
23 #include "Eigen/Core"
24 #include "Eigen/SparseCore"
25 #include "gmock/gmock.h"
26 #include "gtest/gtest.h"
27 #include "ortools/pdlp/quadratic_program.h"
28 #include "ortools/pdlp/sharded_quadratic_program.h"
29 #include "ortools/pdlp/sharder.h"
30 #include "ortools/pdlp/solve_log.pb.h"
31 #include "ortools/pdlp/test_util.h"
32 
33 namespace operations_research::pdlp {
34 namespace {
35 
36 using ::Eigen::VectorXd;
37 using ::testing::ElementsAre;
38 using ::testing::ElementsAreArray;
39 using ::testing::IsNan;
40 
41 TEST(ShardedWeightedAverageTest, SimpleAverage) {
42  Sharder sharder(/*num_elements=*/2, /*num_shards=*/2,
43  /*thread_pool=*/nullptr);
44  Eigen::VectorXd vec1(2), vec2(2);
45  vec1 << 4, 1;
46  vec2 << 1, 7;
47 
48  ShardedWeightedAverage average(&sharder);
49  average.Add(vec1, 1.0);
50  average.Add(vec2, 2.0);
51 
52  ASSERT_TRUE(average.HasNonzeroWeight());
53  EXPECT_EQ(average.NumTerms(), 2);
54 
55  EXPECT_THAT(average.ComputeAverage(), ElementsAre(2.0, 5.0));
56 
57  average.Clear();
58  EXPECT_FALSE(average.HasNonzeroWeight());
59  EXPECT_EQ(average.NumTerms(), 0);
60 }
61 
62 TEST(ShardedWeightedAverageTest, MoveConstruction) {
63  Sharder sharder(/*num_elements=*/2, /*num_shards=*/2,
64  /*thread_pool=*/nullptr);
65  Eigen::VectorXd vec(2);
66  vec << 4, 1;
67 
68  ShardedWeightedAverage average(&sharder);
69  average.Add(vec, 2.0);
70 
71  ShardedWeightedAverage average2(std::move(average));
72  EXPECT_THAT(average2.ComputeAverage(), ElementsAre(4.0, 1.0));
73 }
74 
75 TEST(ShardedWeightedAverageTest, MoveAssignment) {
76  Sharder sharder1(/*num_elements=*/2, /*num_shards=*/2,
77  /*thread_pool=*/nullptr);
78  Sharder sharder2(/*num_elements=*/3, /*num_shards=*/2,
79  /*thread_pool=*/nullptr);
80  Eigen::VectorXd vec1(2), vec2(2);
81  vec1 << 4, 1;
82  vec2 << 0, 3;
83 
84  ShardedWeightedAverage average1(&sharder1);
85  average1.Add(vec1, 2.0);
86 
87  ShardedWeightedAverage average2(&sharder2);
88 
89  average2 = std::move(average1);
90  average2.Add(vec2, 2.0);
91  EXPECT_THAT(average2.ComputeAverage(), ElementsAre(2.0, 2.0));
92 }
93 
94 TEST(ShardedWeightedAverageTest, ZeroAverage) {
95  Sharder sharder(/*num_elements=*/1, /*num_shards=*/1,
96  /*thread_pool=*/nullptr);
97 
98  ShardedWeightedAverage average(&sharder);
99  ASSERT_FALSE(average.HasNonzeroWeight());
100 
101  EXPECT_THAT(average.ComputeAverage(), ElementsAre(0.0));
102 }
103 
104 // This test verifies that if we average an identical vector repeatedly the
105 // average is exactly that vector, with no roundoff.
106 TEST(ShardedWeightedAverageTest, AveragesEqualWithoutRoundoff) {
107  Sharder sharder(/*num_elements=*/4, /*num_shards=*/1,
108  /*thread_pool=*/nullptr);
109  ShardedWeightedAverage average(&sharder);
110  EXPECT_THAT(average.ComputeAverage(), ElementsAre(0, 0, 0, 0));
111  VectorXd data(4);
112  data << 1.0, 1.0 / 3, 3.0 / 7, 3.14159;
113  average.Add(data, 341.45);
114  EXPECT_THAT(average.ComputeAverage(), ElementsAreArray(data));
115  average.Add(data, 1.4134);
116  EXPECT_THAT(average.ComputeAverage(), ElementsAreArray(data));
117  average.Add(data, 7.23);
118  EXPECT_THAT(average.ComputeAverage(), ElementsAreArray(data));
119 }
120 
121 TEST(ShardedWeightedAverageTest, AddsZeroWeight) {
122  Sharder sharder(/*num_elements=*/1, /*num_shards=*/1,
123  /*thread_pool=*/nullptr);
124 
125  ShardedWeightedAverage average(&sharder);
126  ASSERT_FALSE(average.HasNonzeroWeight());
127  VectorXd data(1);
128  data << 1.0;
129  average.Add(data, 0.0);
130  EXPECT_FALSE(average.HasNonzeroWeight());
131  EXPECT_THAT(average.ComputeAverage(), ElementsAre(0.0));
132 }
133 
134 // The combined bounds vector for `TestLp()` is [12, 7, 4, 1].
135 // L_inf norm: 12.0
136 // L_2 norm: sqrt(210.0) ≈ 14.49
137 
138 TEST(ProblemStatsTest, TestLp) {
139  ShardedQuadraticProgram lp(TestLp(), 2, 2);
140  const QuadraticProgramStats stats = ComputeStats(lp);
141 
142  EXPECT_EQ(stats.num_variables(), 4);
143  EXPECT_EQ(stats.num_constraints(), 4);
144  EXPECT_DOUBLE_EQ(stats.constraint_matrix_col_min_l_inf_norm(), 1.0);
145  EXPECT_DOUBLE_EQ(stats.constraint_matrix_row_min_l_inf_norm(), 1.0);
146  EXPECT_EQ(stats.constraint_matrix_num_nonzeros(), 9);
147  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_max(), 4.0);
148  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_min(), 1.0);
149  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_avg(), 14.5 / 9.0);
150  EXPECT_DOUBLE_EQ(stats.constraint_matrix_l2_norm(), std::sqrt(31.25));
151  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_max(), 5.5);
152  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_min(), 1.0);
153  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_avg(), 2.375);
154  EXPECT_DOUBLE_EQ(stats.objective_vector_l2_norm(), std::sqrt(36.25));
155  EXPECT_EQ(stats.objective_matrix_num_nonzeros(), 0);
156  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_max(), 0.0);
157  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_min(), 0.0);
158  EXPECT_THAT(stats.objective_matrix_abs_avg(), IsNan());
159  EXPECT_DOUBLE_EQ(stats.objective_matrix_l2_norm(), 0.0);
160  EXPECT_EQ(stats.variable_bound_gaps_num_finite(), 1);
161  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_max(), 1.0);
162  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_min(), 1.0);
163  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_avg(), 1.0);
164  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_l2_norm(), 1.0);
165  EXPECT_DOUBLE_EQ(stats.combined_bounds_max(), 12.0);
166  EXPECT_DOUBLE_EQ(stats.combined_bounds_min(), 1.0);
167  EXPECT_DOUBLE_EQ(stats.combined_bounds_avg(), 6.0);
168  EXPECT_DOUBLE_EQ(stats.combined_bounds_l2_norm(), std::sqrt(210.0));
169 }
170 
171 TEST(ProblemStatsTest, TinyLp) {
172  ShardedQuadraticProgram lp(TinyLp(), 2, 2);
173  const QuadraticProgramStats stats = ComputeStats(lp);
174 
175  EXPECT_EQ(stats.num_variables(), 4);
176  EXPECT_EQ(stats.num_constraints(), 3);
177  EXPECT_DOUBLE_EQ(stats.constraint_matrix_col_min_l_inf_norm(), 1.0);
178  EXPECT_DOUBLE_EQ(stats.constraint_matrix_row_min_l_inf_norm(), 1.0);
179  EXPECT_EQ(stats.constraint_matrix_num_nonzeros(), 8);
180  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_max(), 2.0);
181  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_min(), 1.0);
182  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_avg(), 1.25);
183  EXPECT_DOUBLE_EQ(stats.constraint_matrix_l2_norm(), std::sqrt(14.0));
184  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_max(), 5.0);
185  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_min(), 1.0);
186  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_avg(), 2.25);
187  EXPECT_DOUBLE_EQ(stats.objective_vector_l2_norm(), std::sqrt(31.0));
188  EXPECT_EQ(stats.objective_matrix_num_nonzeros(), 0);
189  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_max(), 0.0);
190  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_min(), 0.0);
191  EXPECT_THAT(stats.objective_matrix_abs_avg(), IsNan());
192  EXPECT_DOUBLE_EQ(stats.objective_matrix_l2_norm(), 0.0);
193  EXPECT_EQ(stats.variable_bound_gaps_num_finite(), 4);
194  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_max(), 6.0);
195  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_min(), 2.0);
196  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_avg(), 3.75);
197  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_l2_norm(), std::sqrt(65.0));
198  EXPECT_DOUBLE_EQ(stats.combined_bounds_max(), 12.0);
199  EXPECT_DOUBLE_EQ(stats.combined_bounds_min(), 1.0);
200  EXPECT_DOUBLE_EQ(stats.combined_bounds_avg(), 20.0 / 3.0);
201  EXPECT_DOUBLE_EQ(stats.combined_bounds_l2_norm(), std::sqrt(194.0));
202 }
203 
204 TEST(ProblemStatsTest, TestDiagonalQp1) {
205  ShardedQuadraticProgram qp(TestDiagonalQp1(), 2, 2);
206  const QuadraticProgramStats stats = ComputeStats(qp);
207 
208  EXPECT_EQ(stats.num_variables(), 2);
209  EXPECT_EQ(stats.num_constraints(), 1);
210  EXPECT_DOUBLE_EQ(stats.constraint_matrix_col_min_l_inf_norm(), 1.0);
211  EXPECT_DOUBLE_EQ(stats.constraint_matrix_row_min_l_inf_norm(), 1.0);
212  EXPECT_EQ(stats.constraint_matrix_num_nonzeros(), 2);
213  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_max(), 1.0);
214  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_min(), 1.0);
215  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_avg(), 1.0);
216  EXPECT_DOUBLE_EQ(stats.constraint_matrix_l2_norm(), std::sqrt(2.0));
217  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_max(), 1.0);
218  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_min(), 1.0);
219  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_avg(), 1.0);
220  EXPECT_DOUBLE_EQ(stats.objective_vector_l2_norm(), std::sqrt(2.0));
221  EXPECT_EQ(stats.objective_matrix_num_nonzeros(), 2);
222  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_max(), 4.0);
223  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_min(), 1.0);
224  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_avg(), 2.5);
225  EXPECT_DOUBLE_EQ(stats.objective_matrix_l2_norm(), std::sqrt(17.0));
226  EXPECT_EQ(stats.variable_bound_gaps_num_finite(), 2);
227  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_max(), 6.0);
228  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_min(), 1.0);
229  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_avg(), 3.5);
230  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_l2_norm(), std::sqrt(37.0));
231  EXPECT_DOUBLE_EQ(stats.combined_bounds_max(), 1.0);
232  EXPECT_DOUBLE_EQ(stats.combined_bounds_min(), 1.0);
233  EXPECT_DOUBLE_EQ(stats.combined_bounds_avg(), 1.0);
234  EXPECT_DOUBLE_EQ(stats.combined_bounds_l2_norm(), 1.0);
235 }
236 
237 TEST(ProblemStatsTest, ModifiedTestDiagonalQp1) {
238  QuadraticProgram orig_qp = TestDiagonalQp1();
239  // A case where `objective_matrix_num_nonzeros` doesn't match the dimension.
240  orig_qp.objective_matrix->diagonal() << 2.0, 0.0;
241  ShardedQuadraticProgram qp(orig_qp, 2, 2);
242  const QuadraticProgramStats stats = ComputeStats(qp);
243 
244  EXPECT_EQ(stats.objective_matrix_num_nonzeros(), 1);
245  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_max(), 2.0);
246  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_min(), 2.0);
247  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_avg(), 1.0);
248  EXPECT_DOUBLE_EQ(stats.objective_matrix_l2_norm(), 2.0);
249 }
250 
251 // This is like `SmallLp`, except that an `infinite_bound_threshold` of 10
252 // treats the first bound as infinite, leaving [0, 7, 4, 1] as the combined
253 // bounds vector.
254 TEST(ProblemStatsTest, TestLpWithInfiniteConstraintBoundThreshold) {
255  ShardedQuadraticProgram lp(TestLp(), 2, 2);
256  const QuadraticProgramStats stats =
257  ComputeStats(lp, /*infinite_constraint_bound_threshold=*/10);
258 
259  EXPECT_EQ(stats.num_variables(), 4);
260  EXPECT_EQ(stats.num_constraints(), 4);
261  EXPECT_DOUBLE_EQ(stats.constraint_matrix_col_min_l_inf_norm(), 1.0);
262  EXPECT_DOUBLE_EQ(stats.constraint_matrix_row_min_l_inf_norm(), 1.0);
263  EXPECT_EQ(stats.constraint_matrix_num_nonzeros(), 9);
264  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_max(), 4.0);
265  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_min(), 1.0);
266  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_avg(), 14.5 / 9.0);
267  EXPECT_DOUBLE_EQ(stats.constraint_matrix_l2_norm(), std::sqrt(31.25));
268  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_max(), 5.5);
269  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_min(), 1.0);
270  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_avg(), 2.375);
271  EXPECT_DOUBLE_EQ(stats.objective_vector_l2_norm(), std::sqrt(36.25));
272  EXPECT_EQ(stats.objective_matrix_num_nonzeros(), 0);
273  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_max(), 0.0);
274  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_min(), 0.0);
275  EXPECT_THAT(stats.objective_matrix_abs_avg(), IsNan());
276  EXPECT_DOUBLE_EQ(stats.objective_matrix_l2_norm(), 0.0);
277  EXPECT_EQ(stats.variable_bound_gaps_num_finite(), 1);
278  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_max(), 1.0);
279  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_min(), 1.0);
280  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_avg(), 1.0);
281  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_l2_norm(), 1.0);
282  EXPECT_DOUBLE_EQ(stats.combined_bounds_max(), 7.0);
283  EXPECT_DOUBLE_EQ(stats.combined_bounds_min(), 1.0);
284  EXPECT_DOUBLE_EQ(stats.combined_bounds_avg(), 3.0);
285  EXPECT_DOUBLE_EQ(stats.combined_bounds_l2_norm(), std::sqrt(66.0));
286 }
287 
288 TEST(ProblemStatsTest, NoFiniteGaps) {
289  ShardedQuadraticProgram lp(SmallInvalidProblemLp(), 2, 2);
290  const QuadraticProgramStats stats = ComputeStats(lp);
291  // Ensure max/min/avg take their default values when no finite gaps exist.
292  EXPECT_EQ(stats.variable_bound_gaps_num_finite(), 0);
293  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_max(), 0.0);
294  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_min(), 0.0);
295  EXPECT_THAT(stats.variable_bound_gaps_avg(), IsNan());
296  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_l2_norm(), 0.0);
297 }
298 
299 TEST(ProblemStatsTest, LpWithoutConstraints) {
300  ShardedQuadraticProgram lp(LpWithoutConstraints(), 2, 2);
301  const QuadraticProgramStats stats = ComputeStats(lp);
302  // When there are no constraints, max/min absolute values and infinity norms
303  // are assigned 0 by convention. The same is true for the combined bounds.
304  EXPECT_EQ(stats.constraint_matrix_num_nonzeros(), 0);
305  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_max(), 0.0);
306  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_min(), 0.0);
307  EXPECT_THAT(stats.constraint_matrix_abs_avg(), IsNan());
308  EXPECT_DOUBLE_EQ(stats.constraint_matrix_l2_norm(), 0.0);
309  EXPECT_DOUBLE_EQ(stats.constraint_matrix_col_min_l_inf_norm(), 0.0);
310  EXPECT_DOUBLE_EQ(stats.constraint_matrix_row_min_l_inf_norm(), 0.0);
311  EXPECT_DOUBLE_EQ(stats.combined_bounds_max(), 0.0);
312  EXPECT_DOUBLE_EQ(stats.combined_bounds_min(), 0.0);
313  EXPECT_THAT(stats.combined_bounds_avg(), IsNan());
314  EXPECT_DOUBLE_EQ(stats.combined_bounds_l2_norm(), 0.0);
315 }
316 
317 TEST(ProblemStatsTest, EmptyLp) {
318  ShardedQuadraticProgram lp(QuadraticProgram(0, 0), 2, 2);
319  const QuadraticProgramStats stats = ComputeStats(lp);
320  // When `lp` is empty, everything except averages should be 0 and averages
321  // should be NaN.
322  EXPECT_EQ(stats.num_variables(), 0);
323  EXPECT_EQ(stats.num_constraints(), 0);
324  EXPECT_DOUBLE_EQ(stats.constraint_matrix_col_min_l_inf_norm(), 0.0);
325  EXPECT_DOUBLE_EQ(stats.constraint_matrix_row_min_l_inf_norm(), 0.0);
326  EXPECT_EQ(stats.constraint_matrix_num_nonzeros(), 0);
327  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_max(), 0.0);
328  EXPECT_DOUBLE_EQ(stats.constraint_matrix_abs_min(), 0.0);
329  EXPECT_THAT(stats.constraint_matrix_abs_avg(), IsNan());
330  EXPECT_DOUBLE_EQ(stats.constraint_matrix_l2_norm(), 0.0);
331  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_max(), 0.0);
332  EXPECT_DOUBLE_EQ(stats.objective_vector_abs_min(), 0.0);
333  EXPECT_THAT(stats.objective_vector_abs_avg(), IsNan());
334  EXPECT_DOUBLE_EQ(stats.objective_vector_l2_norm(), 0.0);
335  EXPECT_EQ(stats.objective_matrix_num_nonzeros(), 0);
336  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_max(), 0.0);
337  EXPECT_DOUBLE_EQ(stats.objective_matrix_abs_min(), 0.0);
338  EXPECT_THAT(stats.objective_matrix_abs_avg(), IsNan());
339  EXPECT_DOUBLE_EQ(stats.objective_matrix_l2_norm(), 0.0);
340  EXPECT_EQ(stats.variable_bound_gaps_num_finite(), 0);
341  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_max(), 0.0);
342  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_min(), 0.0);
343  EXPECT_THAT(stats.variable_bound_gaps_avg(), IsNan());
344  EXPECT_DOUBLE_EQ(stats.variable_bound_gaps_l2_norm(), 0.0);
345  EXPECT_DOUBLE_EQ(stats.combined_bounds_max(), 0.0);
346  EXPECT_DOUBLE_EQ(stats.combined_bounds_min(), 0.0);
347  EXPECT_THAT(stats.combined_bounds_avg(), IsNan());
348  EXPECT_DOUBLE_EQ(stats.combined_bounds_l2_norm(), 0.0);
349 }
350 
351 // The `TestLp()` matrix is [ 2 1 1 2; 1 0 1 0; 4 0 0 0; 0 0 1.5 -1],
352 // the scaled matrix is [ 0 1 2 -2; 0 0 4 0; 0 0 0 0; 0 0 9 3],
353 // so the row LInf norms are [2 4 0 9] and the column LInf norms are [0 1 9 3].
354 // Rescaling divides the scaling vectors by sqrt(norms).
355 TEST(LInfRuizRescaling, OneIteration) {
356  ShardedQuadraticProgram lp(TestLp(), /*num_threads=*/2, /*num_shards=*/2);
357  VectorXd row_scaling_vec(4), col_scaling_vec(4);
358  row_scaling_vec << 1, 2, 1, 3;
359  col_scaling_vec << 0, 1, 2, -1;
360  LInfRuizRescaling(lp, /*num_iterations=*/1, row_scaling_vec, col_scaling_vec);
361  EXPECT_THAT(row_scaling_vec, ElementsAre(1 / std::sqrt(2), 1.0, 1.0, 1.0));
362  EXPECT_THAT(col_scaling_vec,
363  ElementsAre(0.0, 1.0, 2.0 / 3.0, -1.0 / std::sqrt(3.0)));
364 }
365 
366 // The `TestLp()` matrix is [ 2 1 1 2; 1 0 1 0; 4 0 0 0; 0 0 1.5 -1],
367 // the scaled matrix is [ 0 1 2 -2; 0 0 4 0; 0 0 0 0; 0 0 9 3],
368 // so the row L2 norms are [3 4 0 sqrt(90)] and the column L2 norms are [0 1
369 // sqrt(101) sqrt(13)]. Rescaling divides the scaling vectors by sqrt(norms).
370 TEST(L2RuizRescaling, OneIteration) {
371  ShardedQuadraticProgram lp(TestLp(), /*num_threads=*/2, /*num_shards=*/2);
372  VectorXd row_scaling_vec(4), col_scaling_vec(4);
373  row_scaling_vec << 1, 2, 1, 3;
374  col_scaling_vec << 0, 1, 2, -1;
375  L2NormRescaling(lp, row_scaling_vec, col_scaling_vec);
376  EXPECT_THAT(row_scaling_vec, ElementsAre(1.0 / std::pow(3.0, 0.5), 1.0, 1.0,
377  3.0 / std::pow(90.0, 0.25)));
378  EXPECT_THAT(col_scaling_vec, ElementsAre(0.0, 1.0, 2.0 / std::pow(101, 0.25),
379  -1.0 / std::pow(13.0, 0.25)));
380 }
381 
382 // The `test_lp` matrix is [2 3], so the row L2 norms are [sqrt(13)] and the
383 // column L2 norms are [2 3]. Rescaling divides the scaling vectors by
384 // sqrt(norms).
385 TEST(L2RuizRescaling, OneIterationNonSquare) {
386  QuadraticProgram test_lp(/*num_variables=*/2, /*num_constraints=*/1);
387  std::vector<Eigen::Triplet<double, int64_t>> triplets = {{0, 0, 2.0},
388  {0, 1, 3.0}};
389  test_lp.constraint_matrix.setFromTriplets(triplets.begin(), triplets.end());
390  ShardedQuadraticProgram lp(std::move(test_lp), /*num_threads=*/2,
391  /*num_shards=*/2);
392  VectorXd row_scaling_vec = VectorXd::Ones(1);
393  VectorXd col_scaling_vec = VectorXd::Ones(2);
394  L2NormRescaling(lp, row_scaling_vec, col_scaling_vec);
395  EXPECT_THAT(row_scaling_vec, ElementsAre(1.0 / std::pow(13.0, 0.25)));
396  EXPECT_THAT(col_scaling_vec,
397  ElementsAre(1.0 / std::sqrt(2.0), 1.0 / std::sqrt(3.0)));
398 }
399 
400 // With many iterations of `LInfRuizRescaling`, the scaled matrix should
401 // converge to have col LInf norm 1 and row LInf norm 1.
402 TEST(LInfRuizRescaling, Convergence) {
403  ShardedQuadraticProgram lp(TestLp(), /*num_threads=*/2, /*num_shards=*/2);
404  VectorXd row_scaling_vec(4), col_scaling_vec(4);
405  VectorXd col_norm(4), row_norm(4);
406  row_scaling_vec << 1, 1, 1, 1;
407  col_scaling_vec << 1, 1, 1, 1;
408  LInfRuizRescaling(lp, /*num_iterations=*/20, row_scaling_vec,
409  col_scaling_vec);
410  col_norm = ScaledColLInfNorm(lp.Qp().constraint_matrix, row_scaling_vec,
411  col_scaling_vec, lp.ConstraintMatrixSharder());
412  row_norm = ScaledColLInfNorm(lp.TransposedConstraintMatrix(), col_scaling_vec,
413  row_scaling_vec,
414  lp.TransposedConstraintMatrixSharder());
415  EXPECT_THAT(row_norm, EigenArrayNear<double>({1.0, 1.0, 1.0, 1.0}, 1.0e-4));
416  EXPECT_THAT(col_norm, EigenArrayNear<double>({1.0, 1.0, 1.0, 1.0}, 1.0e-4));
417 }
418 
419 // This applies one round of l_inf and one round of L2 rescaling.
420 // The `TestLp()` matrix is [ 2 1 1 2; 1 0 1 0; 4 0 0 0; 0 0 1.5 -1],
421 // so the row LInf norms are [2 1 4 1.5] and column LInf norms are [4 1 1.5 2].
422 // l_inf divides by sqrt(norms), giving
423 // [0.7071 0.7071 0.5773 1; 0.5 0 0.8165 0; 1 0 0 0; 0 0 1 -0.5773]
424 // which has row L2 norms [1.5275 0.957429 1 1.1547] and col L2 norms
425 // [1.3229 0.7071 1.4142 1.1547]. The resulting scaling vectors are
426 // 1/sqrt((l_inf norms).*(l2 norms)).
427 TEST(ApplyRescaling, ApplyRescalingWorksForTestLp) {
428  ShardedQuadraticProgram lp(TestLp(), /*num_threads=*/2, /*num_shards=*/2);
429  ScalingVectors scaling = ApplyRescaling(
430  RescalingOptions{.l_inf_ruiz_iterations = 1, .l2_norm_rescaling = true},
431  lp);
432  EXPECT_THAT(scaling.row_scaling_vec,
433  EigenArrayNear<double>(
434  {1.0 / sqrt(2.0 * 1.5275), 1.0 / sqrt(1.0 * 0.9574),
435  1.0 / sqrt(4.0 * 1.0), 1.0 / sqrt(1.5 * 1.1547)},
436  1.0e-4));
437  EXPECT_THAT(scaling.col_scaling_vec,
438  EigenArrayNear<double>(
439  {1.0 / sqrt(4.0 * 1.3229), 1.0 / sqrt(1.0 * 0.7071),
440  1.0 / sqrt(1.5 * 1.4142), 1.0 / sqrt(2.0 * 1.1547)},
441  1.0e-4));
442 }
443 
444 TEST(ComputePrimalGradientTest, CorrectForLp) {
445  // The choice of two shards is intentional, to help catch bugs in the sharded
446  // computations.
447  ShardedQuadraticProgram lp(TestLp(), /*num_threads=*/2, /*num_shards=*/2);
448 
449  VectorXd primal_solution(4), dual_solution(4);
450  primal_solution << 0.0, 0.0, 0.0, 3.0;
451  dual_solution << -1.0, 0.0, 1.0, 1.0;
452 
453  const LagrangianPart primal_part = ComputePrimalGradient(
454  lp, primal_solution, lp.TransposedConstraintMatrix() * dual_solution);
455  // Using notation consistent with
456  // https://developers.google.com/optimization/lp/pdlp_math.
457  // c - A^T y
458  EXPECT_THAT(primal_part.gradient,
459  ElementsAre(5.5 - 2.0, -2.0 + 1.0, -1.0 - 0.5, 1.0 + 3.0));
460  // c^T x - y^T Ax.
461  EXPECT_DOUBLE_EQ(primal_part.value, 3.0 + 9.0);
462 }
463 
464 TEST(ComputeDualGradientTest, CorrectForLp) {
465  // The choice of two shards is intentional, to help catch bugs in the sharded
466  // computations.
467  ShardedQuadraticProgram lp(TestLp(), /*num_threads=*/2, /*num_shards=*/2);
468 
469  VectorXd primal_solution(4), dual_solution(4);
470  primal_solution << 0.0, 0.0, 0.0, 3.0;
471  dual_solution << -1.0, 0.0, 1.0, 1.0;
472 
473  const LagrangianPart dual_part = ComputeDualGradient(
474  lp, dual_solution, lp.Qp().constraint_matrix * primal_solution);
475  // Using notation consistent with
476  // https://developers.google.com/optimization/lp/pdlp_math.
477  // active_constraint_right_hand_side - Ax
478  EXPECT_THAT(dual_part.gradient,
479  ElementsAre(12.0 - 6.0, 7.0, -4.0, -1.0 + 3.0));
480  // y^T active_constraint_right_hand_side
481  EXPECT_DOUBLE_EQ(dual_part.value, 12.0 * -1.0 + -4.0 * 1.0 + -1.0 * 1.0);
482 }
483 
484 TEST(ComputeDualGradientTest, CorrectOnTwoSidedConstraints) {
485  QuadraticProgram qp = TestLp();
486  // Makes the constraints all two-sided. The primal solution is feasible in
487  // the first constraint, below the lower bound of the second constraint, and
488  // above the upper bound of the third constraint.
489  qp.constraint_lower_bounds[0] = 4;
490  qp.constraint_lower_bounds[1] = 5;
491  qp.constraint_upper_bounds[2] = -1;
492  ShardedQuadraticProgram sharded_qp(std::move(qp), /*num_threads=*/2,
493  /*num_shards=*/2);
494 
495  VectorXd primal_solution(4), dual_solution(4);
496  primal_solution << 0.0, 0.0, 0.0, 3.0;
497  dual_solution << 0.0, 0.0, 0.0, -1.0;
498 
499  const LagrangianPart dual_part =
500  ComputeDualGradient(sharded_qp, dual_solution,
501  sharded_qp.Qp().constraint_matrix * primal_solution);
502  // Using notation consistent with
503  // https://developers.google.com/optimization/lp/pdlp_math.
504  // active_constraint_right_hand_side - Ax
505  EXPECT_THAT(dual_part.gradient,
506  ElementsAre(0.0, 5.0 - 0.0, -1.0 - 0.0, 1.0 + 3.0));
507  // y^T active_constraint_right_hand_side
508  EXPECT_DOUBLE_EQ(dual_part.value, 1.0 * -1.0);
509 }
510 
511 TEST(HasValidBoundsTest, SmallInvalidLp) {
512  ShardedQuadraticProgram lp(SmallInvalidProblemLp(), /*num_threads=*/2,
513  /*num_shards=*/2);
514 
515  bool is_valid = HasValidBounds(lp);
516  EXPECT_FALSE(is_valid);
517 }
518 
519 TEST(HasValidBoundsTest, SmallValidLp) {
520  ShardedQuadraticProgram lp(SmallPrimalInfeasibleLp(), /*num_threads=*/2,
521  /*num_shards=*/2);
522 
523  bool is_valid = HasValidBounds(lp);
524  EXPECT_TRUE(is_valid);
525 }
526 
527 TEST(ComputePrimalGradientTest, CorrectForQp) {
528  ShardedQuadraticProgram qp(TestDiagonalQp1(), /*num_threads=*/2,
529  /*num_shards=*/2);
530 
531  VectorXd primal_solution(2), dual_solution(1);
532  primal_solution << 1.0, 2.0;
533  dual_solution << -2.0;
534 
535  const LagrangianPart primal_part = ComputePrimalGradient(
536  qp, primal_solution, qp.TransposedConstraintMatrix() * dual_solution);
537 
538  // Using notation consistent with
539  // https://developers.google.com/optimization/lp/pdlp_math.
540  // c - A^T y + Qx
541  EXPECT_THAT(primal_part.gradient,
542  ElementsAre(-1.0 + 2.0 + 4.0, -1.0 + 2.0 + 2.0));
543  // (1/2) x^T Qx + c^T x - y^T Ax.
544  EXPECT_DOUBLE_EQ(primal_part.value, 4.0 - 3.0 + 2.0 * 3.0);
545 }
546 
547 TEST(ComputeDualGradientTest, CorrectForQp) {
548  ShardedQuadraticProgram qp(TestDiagonalQp1(), /*num_threads=*/2,
549  /*num_shards=*/2);
550 
551  VectorXd primal_solution(2), dual_solution(1);
552  primal_solution << 1.0, 2.0;
553  dual_solution << -2.0;
554 
555  const LagrangianPart dual_part = ComputeDualGradient(
556  qp, dual_solution, qp.Qp().constraint_matrix * primal_solution);
557 
558  // Using notation consistent with
559  // https://developers.google.com/optimization/lp/pdlp_math.
560  // active_constraint_right_hand_side - Ax
561  EXPECT_THAT(dual_part.gradient, ElementsAre(1.0 - (1.0 + 2.0)));
562  // y^T active_constraint_right_hand_side
563  EXPECT_DOUBLE_EQ(dual_part.value, -2.0);
564 }
565 
566 TEST(EstimateSingularValuesTest, CorrectForTestLp) {
567  ShardedQuadraticProgram lp(TestLp(), /*num_threads=*/2, /*num_shards=*/2);
568 
569  // The `TestLp()` matrix is [ 2 1 1 2; 1 0 1 0; 4 0 0 0; 0 0 1.5 -1].
570  std::mt19937 random(1);
571  auto result = EstimateMaximumSingularValueOfConstraintMatrix(
572  lp, std::nullopt, std::nullopt,
573  /*desired_relative_error=*/0.01,
574  /*failure_probability=*/0.001, random);
575  EXPECT_NEAR(result.singular_value, 4.76945, 0.01);
576  EXPECT_LT(result.num_iterations, 300);
577 }
578 
579 TEST(EstimateSingularValuesTest, CorrectForTestLpWithActivePrimalSubspace) {
580  ShardedQuadraticProgram lp(TestLp(), /*num_threads=*/2, /*num_shards=*/2);
581 
582  VectorXd primal_solution(4);
583  // Chosen so x_1 is at its bound, and all other variables are not at bounds.
584  primal_solution << 0.0, -2.0, 0.0, 3.0;
585  // The `TestLp()` matrix is [ 2 1 1 2; 1 0 1 0; 4 0 0 0; 0 0 1.5 -1],
586  // so the projected matrix is [ 2 1 2; 1 1 0; 4 0 0; 0 1.5 -1].
587  std::mt19937 random(1);
588  auto result = EstimateMaximumSingularValueOfConstraintMatrix(
589  lp, primal_solution, std::nullopt, /*desired_relative_error=*/0.01,
590  /*failure_probability=*/0.001, random);
591  EXPECT_NEAR(result.singular_value, 4.73818, 0.01);
592  EXPECT_LT(result.num_iterations, 300);
593 }
594 
595 TEST(EstimateSingularValuesTest, CorrectForTestLpWithActiveDualSubspace) {
596  ShardedQuadraticProgram lp(TestLp(), /*num_threads=*/2, /*num_shards=*/2);
597 
598  VectorXd dual_solution(4);
599  // Chosen so the second dual is at its bound, and all other duals are not at
600  // bounds.
601  dual_solution << 1.0, 0.0, 1.0, 3.0;
602  // The `TestLp()` matrix is [ 2 1 1 2; 1 0 1 0; 4 0 0 0; 0 0 1.5 -1],
603  // so the projected matrix is [ 2 1 1 2; 4 0 0 0; 0 0 1.5 -1].
604  std::mt19937 random(1);
605  auto result = EstimateMaximumSingularValueOfConstraintMatrix(
606  lp, std::nullopt, dual_solution, /*desired_relative_error=*/0.01,
607  /*failure_probability=*/0.001, random);
608  EXPECT_NEAR(result.singular_value, 4.64203, 0.01);
609  EXPECT_LT(result.num_iterations, 300);
610 }
611 
612 TEST(EstimateSingularValuesTest, CorrectForTestLpWithBothActiveSubspaces) {
613  ShardedQuadraticProgram lp(TestLp(), /*num_threads=*/2, /*num_shards=*/2);
614 
615  VectorXd primal_solution(4), dual_solution(4);
616  // Chosen so x_1 is at its bound, and all other variables are not at bounds.
617  primal_solution << 0.0, -2.0, 0.0, 3.0;
618  // Chosen so the second dual is at its bound, and all other duals are not at
619  // bounds.
620  dual_solution << 1.0, 0.0, 1.0, 3.0;
621  // The `TestLp()` matrix is [ 2 1 1 2; 1 0 1 0; 4 0 0 0; 0 0 1.5 -1],
622  // so the projected matrix is [ 2 1 2; 4 0 0; 0 1.5 -1].
623  std::mt19937 random(1);
624  auto result = EstimateMaximumSingularValueOfConstraintMatrix(
625  lp, primal_solution, dual_solution, /*desired_relative_error=*/0.01,
626  /*failure_probability=*/0.001, random);
627  EXPECT_NEAR(result.singular_value, 4.60829, 0.01);
628  EXPECT_LT(result.num_iterations, 300);
629 }
630 
631 TEST(EstimateSingularValuesTest, CorrectForDiagonalLp) {
632  QuadraticProgram diagonal_lp = TestLp();
633  std::vector<Eigen::Triplet<double, int64_t>> triplets = {
634  {0, 0, 2}, {1, 1, 1}, {2, 2, -3}, {3, 3, -1}};
635  diagonal_lp.constraint_matrix.setFromTriplets(triplets.begin(),
636  triplets.end());
637  ShardedQuadraticProgram lp(diagonal_lp, /*num_threads=*/2, /*num_shards=*/2);
638 
639  // The `diagonal_lp` matrix is [ 2 0 0 0; 0 1 0 0; 0 0 -3 0; 0 0 0 -1].
640  std::mt19937 random(1);
641  auto result = EstimateMaximumSingularValueOfConstraintMatrix(
642  lp, std::nullopt, std::nullopt,
643  /*desired_relative_error=*/0.01,
644  /*failure_probability=*/0.001, random);
645  EXPECT_NEAR(result.singular_value, 3, 0.0001);
646  EXPECT_LT(result.num_iterations, 300);
647 }
648 
649 TEST(ProjectToPrimalVariableBoundsTest, TestLp) {
650  ShardedQuadraticProgram qp(TestLp(), /*num_threads=*/2,
651  /*num_shards=*/2);
652  VectorXd primal(4);
653  primal << -3, -3, 5, 5;
654  ProjectToPrimalVariableBounds(qp, primal);
655  EXPECT_THAT(primal, ElementsAre(-3, -2, 5, 3.5));
656 }
657 
658 TEST(ProjectToDualVariableBoundsTest, TestLp) {
659  ShardedQuadraticProgram qp(TestLp(), /*num_threads=*/2,
660  /*num_shards=*/2);
661  VectorXd dual(4);
662  dual << 1, 1, -1, -1;
663  ProjectToDualVariableBounds(qp, dual);
664  EXPECT_THAT(dual, ElementsAre(1, 0, 0, -1));
665 }
666 
667 } // namespace
668 } // namespace operations_research::pdlp
operations_research::pdlp
Definition: iteration_stats.cc:35
operations_research::pdlp::ComputeDualGradient
LagrangianPart ComputeDualGradient(const ShardedQuadraticProgram &sharded_qp, const VectorXd &dual_solution, const VectorXd &primal_product)
Definition: sharded_optimization_utils.cc:497
operations_research::pdlp::ComputePrimalGradient
LagrangianPart ComputePrimalGradient(const ShardedQuadraticProgram &sharded_qp, const VectorXd &primal_solution, const VectorXd &dual_product)
Definition: sharded_optimization_utils.cc:441
operations_research::pdlp::LInfRuizRescaling
void LInfRuizRescaling(const ShardedQuadraticProgram &sharded_qp, const int num_iterations, VectorXd &row_scaling_vec, VectorXd &col_scaling_vec)
Definition: sharded_optimization_utils.cc:404
operations_research::pdlp::EstimateMaximumSingularValueOfConstraintMatrix
SingularValueAndIterations EstimateMaximumSingularValueOfConstraintMatrix(const ShardedQuadraticProgram &sharded_qp, const std::optional< VectorXd > &primal_solution, const std::optional< VectorXd > &dual_solution, const double desired_relative_error, const double failure_probability, std::mt19937 &mt_generator)
Definition: sharded_optimization_utils.cc:671
operations_research::pdlp::ScaledColLInfNorm
VectorXd ScaledColLInfNorm(const Eigen::SparseMatrix< double, Eigen::ColMajor, int64_t > &matrix, const VectorXd &row_scaling_vec, const VectorXd &col_scaling_vec, const Sharder &sharder)
Definition: sharder.cc:286
operations_research::pdlp::HasValidBounds
bool HasValidBounds(const QuadraticProgram &qp)
Definition: quadratic_program.cc:82
operations_research::pdlp::TinyLp
QuadraticProgram TinyLp()
Definition: test_util.cc:67
operations_research::pdlp::ProjectToDualVariableBounds
void ProjectToDualVariableBounds(const ShardedQuadraticProgram &sharded_qp, VectorXd &dual)
Definition: sharded_optimization_utils.cc:727
operations_research::pdlp::LpWithoutConstraints
QuadraticProgram LpWithoutConstraints()
Definition: test_util.cc:262
operations_research::pdlp::SmallInvalidProblemLp
QuadraticProgram SmallInvalidProblemLp()
Definition: test_util.cc:191
operations_research::pdlp::L2NormRescaling
void L2NormRescaling(const ShardedQuadraticProgram &sharded_qp, VectorXd &row_scaling_vec, VectorXd &col_scaling_vec)
Definition: sharded_optimization_utils.cc:411
operations_research::pdlp::ApplyRescaling
ScalingVectors ApplyRescaling(const RescalingOptions &rescaling_options, ShardedQuadraticProgram &sharded_qp)
Definition: sharded_optimization_utils.cc:418
operations_research::pdlp::ComputeStats
QuadraticProgramStats ComputeStats(const ShardedQuadraticProgram &qp, const double infinite_constraint_bound_threshold)
Definition: sharded_optimization_utils.cc:274
operations_research::pdlp::ProjectToPrimalVariableBounds
void ProjectToPrimalVariableBounds(const ShardedQuadraticProgram &sharded_qp, VectorXd &primal)
Definition: sharded_optimization_utils.cc:716
operations_research::pdlp::TestLp
QuadraticProgram TestLp()
Definition: test_util.cc:33
operations_research::pdlp::SmallPrimalInfeasibleLp
QuadraticProgram SmallPrimalInfeasibleLp()
Definition: test_util.cc:217
operations_research::pdlp::TestDiagonalQp1
QuadraticProgram TestDiagonalQp1()
Definition: test_util.cc:143
operations_research::TEST
TEST(LinearAssignmentTest, NullMatrix)
Definition: hungarian_test.cc:78
quadratic_program.h
average
double average
Definition: sharded_optimization_utils.cc:104
sharded_optimization_utils.h
sharded_quadratic_program.h
test_util.h