Line data Source code
1 : /**
2 : * @file test_QuadricProblem.cpp
3 : *
4 : * @brief Tests for the QuqricProblem class
5 : *
6 : * @author Nikola Dinev
7 : */
8 :
9 : #include "doctest/doctest.h"
10 :
11 : #include "QuadricProblem.h"
12 : #include "Identity.h"
13 : #include "Scaling.h"
14 : #include "L2NormPow2.h"
15 : #include "WeightedL2NormPow2.h"
16 : #include "Huber.h"
17 : #include "Logger.h"
18 : #include "VolumeDescriptor.h"
19 : #include "testHelpers.h"
20 : #include "TypeCasts.hpp"
21 :
22 : #include <array>
23 :
24 : using namespace elsa;
25 : using namespace doctest;
26 :
27 : template <template <typename> typename T, typename data_t>
28 : constexpr data_t return_data_t(const T<data_t>&);
29 :
30 : TEST_SUITE_BEGIN("problems");
31 :
32 : TEST_CASE_TEMPLATE("QuadricProblem: Construction with a Quadric functional", TestType, float,
33 : double)
34 8 : {
35 8 : using data_t = TestType;
36 :
37 : // eliminate the timing info from console for the tests
38 8 : Logger::setLevel(Logger::LogLevel::WARN);
39 :
40 8 : GIVEN("a Quadric functional")
41 8 : {
42 8 : IndexVector_t numCoeff(3);
43 8 : numCoeff << 13, 11, 7;
44 8 : VolumeDescriptor dd{numCoeff};
45 :
46 8 : auto scaleFactor = static_cast<data_t>(3.0);
47 8 : Scaling<data_t> scalingOp{dd, scaleFactor};
48 :
49 8 : Eigen::Matrix<data_t, -1, 1> randomData{dd.getNumberOfCoefficients()};
50 8 : randomData.setRandom();
51 8 : DataContainer<data_t> dc{dd, randomData};
52 :
53 8 : Quadric<data_t> quad{scalingOp, dc};
54 :
55 8 : WHEN("constructing a QuadricProblem without x0 from it")
56 8 : {
57 4 : QuadricProblem<data_t> prob{quad};
58 :
59 4 : THEN("the clone works correctly")
60 4 : {
61 2 : auto probClone = prob.clone();
62 :
63 2 : REQUIRE_NE(probClone.get(), &prob);
64 2 : REQUIRE_EQ(*probClone, prob);
65 2 : }
66 :
67 4 : THEN("the problem behaves as expected")
68 4 : {
69 2 : DataContainer<data_t> dcZero(dd);
70 2 : dcZero = 0;
71 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), dcZero));
72 :
73 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(), 0));
74 2 : REQUIRE_UNARY(checkApproxEq(prob.getGradient(), as<data_t>(-1.0) * dc));
75 :
76 2 : auto hessian = prob.getHessian();
77 2 : REQUIRE_EQ(hessian, leaf(scalingOp));
78 2 : }
79 4 : }
80 :
81 8 : WHEN("constructing a QuadricProblem with x0 from it")
82 8 : {
83 4 : Eigen::Matrix<data_t, -1, 1> x0Vec{dd.getNumberOfCoefficients()};
84 4 : x0Vec.setRandom();
85 4 : DataContainer<data_t> x0{dd, x0Vec};
86 :
87 4 : QuadricProblem<data_t> prob{quad, x0};
88 :
89 4 : THEN("the clone works correctly")
90 4 : {
91 2 : auto probClone = prob.clone();
92 :
93 2 : REQUIRE_NE(probClone.get(), &prob);
94 2 : REQUIRE_EQ(*probClone, prob);
95 2 : }
96 :
97 4 : THEN("the problem behaves as expected")
98 4 : {
99 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
100 :
101 2 : REQUIRE_UNARY(checkApproxEq(
102 2 : prob.evaluate(),
103 2 : static_cast<data_t>(0.5 * scaleFactor) * x0.squaredL2Norm() - x0.dot(dc)));
104 :
105 2 : auto gradient = prob.getGradient();
106 2 : DataContainer gradientDirect = scaleFactor * x0 - dc;
107 2004 : for (index_t i = 0; i < gradient.getSize(); ++i)
108 2 : REQUIRE_UNARY(checkApproxEq(gradient[i], gradientDirect[i]));
109 :
110 2 : auto hessian = prob.getHessian();
111 2 : REQUIRE_EQ(hessian, leaf(scalingOp));
112 2 : }
113 4 : }
114 8 : }
115 8 : }
116 :
117 : TEST_CASE_TEMPLATE("QuadricProblem: with a Quadric functional with spd operator", TestType, float,
118 : double)
119 8 : {
120 8 : using data_t = TestType;
121 :
122 : // eliminate the timing info from console for the tests
123 8 : Logger::setLevel(Logger::LogLevel::WARN);
124 :
125 8 : GIVEN("a spd operator")
126 8 : {
127 8 : IndexVector_t numCoeff(3);
128 8 : numCoeff << 13, 11, 7;
129 8 : VolumeDescriptor dd{numCoeff};
130 :
131 8 : data_t scaleFactor = 3.0;
132 8 : Scaling<data_t> scalingOp{dd, scaleFactor};
133 :
134 8 : Eigen::Matrix<data_t, -1, 1> randomData{dd.getNumberOfCoefficients()};
135 8 : randomData.setRandom();
136 8 : DataContainer<data_t> dc{dd, randomData};
137 :
138 8 : WHEN("constructing a QuadricProblem without x0 from it")
139 8 : {
140 4 : QuadricProblem<data_t> prob{scalingOp, dc, true};
141 :
142 4 : THEN("the clone works correctly")
143 4 : {
144 2 : auto probClone = prob.clone();
145 :
146 2 : REQUIRE_NE(probClone.get(), &prob);
147 2 : REQUIRE_EQ(*probClone, prob);
148 2 : }
149 :
150 4 : THEN("the problem behaves as expected")
151 4 : {
152 2 : DataContainer<data_t> dcZero(dd);
153 2 : dcZero = 0;
154 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), dcZero));
155 :
156 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(), 0));
157 2 : REQUIRE_UNARY(isApprox(prob.getGradient(), static_cast<data_t>(-1.0) * dc));
158 :
159 2 : auto hessian = prob.getHessian();
160 2 : REQUIRE_EQ(hessian, leaf(scalingOp));
161 2 : }
162 4 : }
163 :
164 8 : WHEN("constructing a QuadricProblem with x0 from it")
165 8 : {
166 4 : Eigen::Matrix<data_t, -1, 1> x0Vec{dd.getNumberOfCoefficients()};
167 4 : x0Vec.setRandom();
168 4 : DataContainer<data_t> x0{dd, x0Vec};
169 :
170 4 : QuadricProblem<data_t> prob{scalingOp, dc, x0, true};
171 :
172 4 : THEN("the clone works correctly")
173 4 : {
174 2 : auto probClone = prob.clone();
175 :
176 2 : REQUIRE_NE(probClone.get(), &prob);
177 2 : REQUIRE_EQ(*probClone, prob);
178 2 : }
179 :
180 4 : THEN("the problem behaves as expected")
181 4 : {
182 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
183 :
184 2 : REQUIRE_UNARY(
185 2 : checkApproxEq(prob.evaluate(),
186 2 : as<data_t>(0.5 * scaleFactor) * x0.squaredL2Norm() - x0.dot(dc)));
187 :
188 2 : auto gradient = prob.getGradient();
189 2 : DataContainer gradientDirect = scaleFactor * x0 - dc;
190 2004 : for (index_t i = 0; i < gradient.getSize(); ++i)
191 2 : REQUIRE_UNARY(checkApproxEq(gradient[i], gradientDirect[i]));
192 :
193 2 : auto hessian = prob.getHessian();
194 2 : REQUIRE_EQ(hessian, leaf(scalingOp));
195 2 : }
196 4 : }
197 8 : }
198 8 : }
199 :
200 : TEST_CASE_TEMPLATE("QuadricProblem: with a Quadric functional with non-spd operator", TestType,
201 : float, double)
202 8 : {
203 8 : using data_t = TestType;
204 :
205 : // eliminate the timing info from console for the tests
206 8 : Logger::setLevel(Logger::LogLevel::WARN);
207 :
208 8 : GIVEN("a non-spd operator")
209 8 : {
210 8 : IndexVector_t numCoeff(3);
211 8 : numCoeff << 13, 11, 7;
212 8 : VolumeDescriptor dd{numCoeff};
213 :
214 8 : data_t scaleFactor = -3.0;
215 8 : Scaling<data_t> scalingOp{dd, scaleFactor};
216 :
217 8 : Eigen::Matrix<data_t, -1, 1> randomData{dd.getNumberOfCoefficients()};
218 8 : randomData.setRandom();
219 8 : DataContainer<data_t> dc{dd, randomData};
220 :
221 8 : WHEN("Constructing a QuadricProblem without x0 from it")
222 8 : {
223 4 : QuadricProblem<data_t> prob{scalingOp, dc, false};
224 :
225 4 : THEN("the clone works correctly")
226 4 : {
227 2 : auto probClone = prob.clone();
228 :
229 2 : REQUIRE_NE(probClone.get(), &prob);
230 2 : REQUIRE_EQ(*probClone, prob);
231 2 : }
232 :
233 4 : THEN("the problem behaves as expected")
234 4 : {
235 2 : DataContainer<data_t> dcZero(dd);
236 2 : dcZero = 0;
237 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), dcZero));
238 :
239 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(), 0));
240 2 : REQUIRE_UNARY(isApprox(prob.getGradient(), as<data_t>(-scaleFactor) * dc));
241 :
242 2 : auto hessian = prob.getHessian();
243 2 : REQUIRE_EQ(hessian, leaf(adjoint(scalingOp) * scalingOp));
244 2 : }
245 4 : }
246 :
247 8 : WHEN("Constructing a QuadricProblem with x0 from it")
248 8 : {
249 4 : Eigen::Matrix<data_t, -1, 1> x0Vec{dd.getNumberOfCoefficients()};
250 4 : x0Vec.setRandom();
251 4 : DataContainer<data_t> x0{dd, x0Vec};
252 :
253 4 : QuadricProblem<data_t> prob{scalingOp, dc, x0, false};
254 :
255 4 : THEN("the clone works correctly")
256 4 : {
257 2 : auto probClone = prob.clone();
258 :
259 2 : REQUIRE_NE(probClone.get(), &prob);
260 2 : REQUIRE_EQ(*probClone, prob);
261 2 : }
262 :
263 4 : THEN("the problem behaves as expected")
264 4 : {
265 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
266 :
267 2 : REQUIRE_UNARY(
268 2 : checkApproxEq(prob.evaluate(),
269 2 : as<data_t>(0.5 * scaleFactor * scaleFactor) * x0.squaredL2Norm()
270 2 : - scaleFactor * x0.dot(dc)));
271 :
272 2 : auto gradient = prob.getGradient();
273 2 : DataContainer gradientDirect = scaleFactor * (scaleFactor * x0) - scaleFactor * dc;
274 2004 : for (index_t i = 0; i < gradient.getSize(); ++i)
275 2 : REQUIRE_UNARY(checkApproxEq(gradient[i], gradientDirect[i]));
276 :
277 2 : auto hessian = prob.getHessian();
278 2 : REQUIRE_EQ(hessian, leaf(adjoint(scalingOp) * scalingOp));
279 2 : }
280 4 : }
281 8 : }
282 8 : }
283 :
284 : TEST_CASE_TEMPLATE("QuadricProblem: with a different optimization problems", TestType, float,
285 : double)
286 26 : {
287 26 : using data_t = TestType;
288 :
289 : // eliminate the timing info from console for the tests
290 26 : Logger::setLevel(Logger::LogLevel::WARN);
291 :
292 26 : GIVEN("an optimization problem")
293 26 : {
294 26 : IndexVector_t numCoeff(3);
295 26 : numCoeff << 13, 11, 7;
296 26 : VolumeDescriptor dd{numCoeff};
297 :
298 26 : data_t scaleFactor = -3.0;
299 26 : Scaling<data_t> scalingOp{dd, scaleFactor};
300 :
301 26 : Eigen::Matrix<data_t, -1, 1> randomData{dd.getNumberOfCoefficients()};
302 26 : randomData.setRandom();
303 26 : DataContainer<data_t> dc{dd, randomData};
304 :
305 26 : data_t weightFactor = 2.0;
306 26 : Scaling<data_t> weightingOp{dd, weightFactor};
307 :
308 26 : WHEN("trying to convert a problem with a non-quadric non-wls data term")
309 26 : {
310 2 : Problem prob{Huber<data_t>{dd}};
311 2 : THEN("an exception is thrown") { REQUIRE_THROWS(QuadricProblem{prob}); }
312 2 : }
313 :
314 26 : WHEN("trying to convert a problem with a non-quadric non-wls regularization term")
315 26 : {
316 2 : Problem prob{L2NormPow2<data_t>{dd},
317 2 : RegularizationTerm{static_cast<data_t>(0.5), Huber<data_t>{dd}}};
318 2 : THEN("an exception is thrown") { REQUIRE_THROWS(QuadricProblem{prob}); }
319 2 : }
320 :
321 26 : WHEN("converting an optimization problem that has only a quadric data term")
322 26 : {
323 4 : randomData.setRandom();
324 4 : DataContainer<data_t> x0{dd, randomData};
325 :
326 4 : Problem<data_t> initialProb{Quadric<data_t>{weightingOp, dc}, x0};
327 :
328 4 : QuadricProblem<data_t> prob{initialProb};
329 4 : THEN("the clone works correctly")
330 4 : {
331 2 : auto probClone = prob.clone();
332 :
333 2 : REQUIRE_NE(probClone.get(), &prob);
334 2 : REQUIRE_EQ(*probClone, prob);
335 2 : }
336 :
337 4 : THEN("the problem behaves as expected")
338 4 : {
339 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
340 :
341 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
342 2 : as<data_t>(0.5 * weightFactor) * x0.squaredL2Norm()
343 2 : - x0.dot(dc)));
344 2 : REQUIRE_UNARY(isApprox(prob.getGradient(), (weightFactor * x0) - dc));
345 :
346 2 : auto hessian = prob.getHessian();
347 2 : REQUIRE_EQ(hessian, leaf(weightingOp));
348 2 : }
349 4 : }
350 :
351 26 : WHEN("calling the conversion constructor on a QadricProblem")
352 26 : {
353 6 : randomData.setRandom();
354 6 : DataContainer<data_t> x0{dd, randomData};
355 :
356 6 : QuadricProblem<data_t> initialProb{Quadric<data_t>{weightingOp, dc}, x0};
357 :
358 6 : QuadricProblem<data_t> prob{initialProb};
359 6 : THEN("conversion yields the same result as cloning")
360 6 : {
361 2 : REQUIRE_EQ(*initialProb.clone(), prob);
362 2 : }
363 :
364 6 : THEN("the clone works correctly")
365 6 : {
366 2 : auto probClone = prob.clone();
367 :
368 2 : REQUIRE_NE(probClone.get(), &prob);
369 2 : REQUIRE_EQ(*probClone, prob);
370 2 : }
371 :
372 6 : THEN("the problem behaves as expected")
373 6 : {
374 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
375 :
376 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
377 2 : as<data_t>(0.5 * weightFactor) * x0.squaredL2Norm()
378 2 : - x0.dot(dc)));
379 2 : REQUIRE_UNARY(isApprox(prob.getGradient(), (weightFactor * x0) - dc));
380 :
381 2 : auto hessian = prob.getHessian();
382 2 : REQUIRE_EQ(hessian, leaf(weightingOp));
383 2 : }
384 6 : }
385 :
386 26 : randomData.setRandom();
387 26 : DataContainer<data_t> x0{dd, randomData};
388 :
389 26 : Scaling A{dd, static_cast<data_t>(2.0)};
390 26 : randomData.setRandom();
391 26 : DataContainer<data_t> b{dd, randomData};
392 :
393 26 : std::vector<std::unique_ptr<Functional<data_t>>> dataTerms;
394 26 : dataTerms.push_back(std::make_unique<Quadric<data_t>>(dd));
395 26 : dataTerms.push_back(std::make_unique<Quadric<data_t>>(A));
396 26 : dataTerms.push_back(std::make_unique<Quadric<data_t>>(b));
397 26 : dataTerms.push_back(std::make_unique<Quadric<data_t>>(A, b));
398 :
399 26 : Scaling L{dd, static_cast<data_t>(0.5)};
400 26 : randomData.setRandom();
401 26 : DataContainer c{dd, randomData};
402 :
403 26 : std::vector<std::unique_ptr<Quadric<data_t>>> regFunc;
404 26 : regFunc.push_back(std::make_unique<Quadric<data_t>>(dd));
405 26 : regFunc.push_back(std::make_unique<Quadric<data_t>>(A));
406 26 : regFunc.push_back(std::make_unique<Quadric<data_t>>(b));
407 26 : regFunc.push_back(std::make_unique<Quadric<data_t>>(A, b));
408 :
409 26 : std::array descriptions = {"has no operator and no vector", "has an operator but no vector",
410 26 : "has no operator, but has a vector",
411 26 : "has an operator and a vector"};
412 :
413 130 : for (std::size_t i = 0; i < dataTerms.size(); i++) {
414 520 : for (std::size_t j = 0; j < regFunc.size(); j++) {
415 416 : WHEN("Calling the conversion constructor with an OptimizationProblem with only "
416 416 : "quadric terms")
417 416 : {
418 2 : INFO("The data term: ", descriptions[i]);
419 2 : INFO("The regularization term: ", descriptions[j]);
420 :
421 2 : RegularizationTerm reg{static_cast<data_t>(0.25), *regFunc[j]};
422 2 : Problem prob{*dataTerms[i], reg, x0};
423 :
424 2 : QuadricProblem converted{prob};
425 :
426 2 : THEN("the problem can be converted and all operations yield the same result as "
427 2 : "for the initial problem")
428 2 : {
429 :
430 2 : REQUIRE_UNARY(checkApproxEq(converted.evaluate(), prob.evaluate()));
431 :
432 2 : DataContainer<data_t> gradDiff =
433 2 : prob.getGradient() - converted.getGradient();
434 2 : REQUIRE_UNARY(checkApproxEq(gradDiff.squaredL2Norm(), 0));
435 :
436 2 : REQUIRE_UNARY(isApprox(prob.getHessian().apply(x0),
437 2 : converted.getHessian().apply(x0)));
438 2 : }
439 2 : }
440 416 : }
441 104 : }
442 :
443 26 : dataTerms.clear();
444 26 : dataTerms.push_back(std::move(std::make_unique<L2NormPow2<data_t>>(dd)));
445 26 : dataTerms.push_back(
446 26 : std::move(std::make_unique<L2NormPow2<data_t>>(LinearResidual<data_t>{scalingOp})));
447 26 : dataTerms.push_back(
448 26 : std::move(std::make_unique<L2NormPow2<data_t>>(LinearResidual<data_t>{dc})));
449 26 : dataTerms.push_back(
450 26 : std::move(std::make_unique<L2NormPow2<data_t>>(LinearResidual<data_t>{scalingOp, dc})));
451 26 : dataTerms.push_back(std::move(std::make_unique<WeightedL2NormPow2<data_t>>(weightingOp)));
452 26 : dataTerms.push_back(std::move(std::make_unique<WeightedL2NormPow2<data_t>>(
453 26 : LinearResidual<data_t>{scalingOp}, weightingOp)));
454 26 : dataTerms.push_back(std::move(
455 26 : std::make_unique<WeightedL2NormPow2<data_t>>(LinearResidual<data_t>{dc}, weightingOp)));
456 26 : dataTerms.push_back(std::move(std::make_unique<WeightedL2NormPow2<data_t>>(
457 26 : LinearResidual<data_t>{scalingOp, dc}, weightingOp)));
458 :
459 208 : for (const auto& dataTerm : dataTerms) {
460 208 : const auto& res = static_cast<const LinearResidual<data_t>&>(dataTerm->getResidual());
461 208 : const auto isWeighted = is<WeightedL2NormPow2<data_t>>(dataTerm.get());
462 208 : std::string probType = isWeighted ? "wls problem" : "ls problem";
463 :
464 208 : std::string desc =
465 208 : probType
466 208 : + (std::string(res.hasOperator() ? " with an operator" : " with no operator")
467 208 : + ", ")
468 208 : + (res.hasDataVector() ? "a data vector" : "no data vector");
469 208 : WHEN("Converting a quadric problem with no x0")
470 208 : {
471 2 : INFO("Converting a ", desc, "and no x0 to a quadric problem");
472 : // use OptimizationProblem istead of WLSProblem as it allows for more general
473 : // formulations
474 2 : Problem<data_t> initialProb{*dataTerm};
475 :
476 2 : QuadricProblem<data_t> prob{initialProb};
477 :
478 2 : THEN("the clone works correctly")
479 2 : {
480 2 : auto probClone = prob.clone();
481 :
482 2 : REQUIRE_NE(probClone.get(), &prob);
483 2 : REQUIRE_EQ(*probClone, prob);
484 :
485 2 : AND_THEN("the problem behaves as expected")
486 2 : {
487 2 : DataContainer<data_t> zero{dd};
488 2 : zero = 0;
489 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), zero));
490 :
491 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(), 0));
492 :
493 2 : DataContainer<data_t> grad =
494 2 : res.hasDataVector() ? static_cast<data_t>(-1.0) * dc : zero;
495 2 : if (res.hasDataVector() && res.hasOperator())
496 0 : grad *= scaleFactor;
497 2 : if (res.hasDataVector() && isWeighted)
498 0 : grad *= weightFactor;
499 2 : REQUIRE_UNARY(isApprox(prob.getGradient(), grad));
500 :
501 2 : if (isWeighted) {
502 0 : if (res.hasOperator()) {
503 0 : REQUIRE_EQ(prob.getHessian(),
504 0 : leaf(adjoint(scalingOp) * weightingOp * scalingOp));
505 0 : } else {
506 0 : REQUIRE_EQ(prob.getHessian(), leaf(weightingOp));
507 0 : }
508 2 : } else {
509 2 : if (res.hasOperator()) {
510 0 : REQUIRE_EQ(prob.getHessian(), leaf(adjoint(scalingOp) * scalingOp));
511 2 : } else {
512 2 : REQUIRE_EQ(prob.getHessian(), leaf(Identity<data_t>{dd}));
513 2 : }
514 2 : }
515 2 : }
516 2 : }
517 2 : }
518 :
519 208 : WHEN("Converting a quadric problem with x0")
520 208 : {
521 2 : INFO("Converting a ", desc, "and x0 to a quadric problem");
522 2 : randomData.setRandom();
523 2 : DataContainer<data_t> x0{dd, randomData};
524 :
525 : // use OptimizationProblem istead of WLSProblem as it allows for more general
526 : // formulations
527 2 : Problem<data_t> initialProb{*dataTerm, x0};
528 :
529 2 : QuadricProblem<data_t> prob{initialProb};
530 :
531 2 : THEN("the clone works correctly")
532 2 : {
533 2 : auto probClone = prob.clone();
534 :
535 2 : REQUIRE_NE(probClone.get(), &prob);
536 2 : REQUIRE_EQ(*probClone, prob);
537 :
538 2 : AND_THEN("the problem behaves as expected")
539 2 : {
540 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
541 :
542 2 : if (isWeighted) {
543 0 : if (res.hasOperator()) {
544 0 : if (res.hasDataVector()) {
545 0 : REQUIRE_UNARY(
546 0 : checkApproxEq(prob.evaluate(),
547 0 : (as<data_t>(0.5) * scaleFactor * scaleFactor
548 0 : * weightFactor * x0.squaredL2Norm()
549 0 : - scaleFactor * weightFactor * x0.dot(dc))));
550 0 : auto gradient = prob.getGradient();
551 0 : DataContainer gradientDirect =
552 0 : scaleFactor * (weightFactor * (scaleFactor * x0))
553 0 : - scaleFactor * (weightFactor * dc);
554 0 : for (index_t i = 0; i < gradient.getSize(); ++i)
555 0 : REQUIRE_UNARY(
556 0 : checkApproxEq(gradient[i], gradientDirect[i]));
557 0 : } else {
558 0 : REQUIRE_UNARY(checkApproxEq(
559 0 : prob.evaluate(), as<data_t>(0.5) * scaleFactor * scaleFactor
560 0 : * weightFactor * x0.squaredL2Norm()));
561 :
562 0 : auto gradient = prob.getGradient();
563 0 : DataContainer gradientDirect =
564 0 : scaleFactor * (weightFactor * (scaleFactor * x0));
565 0 : for (index_t i = 0; i < gradient.getSize(); ++i)
566 0 : REQUIRE_UNARY(
567 0 : checkApproxEq(gradient[i], gradientDirect[i]));
568 0 : }
569 0 : REQUIRE_EQ(prob.getHessian(),
570 0 : leaf(adjoint(scalingOp) * weightingOp * scalingOp));
571 0 : } else {
572 0 : if (res.hasDataVector()) {
573 0 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
574 0 : as<data_t>(0.5) * weightFactor
575 0 : * x0.squaredL2Norm()
576 0 : - weightFactor * x0.dot(dc)));
577 0 : REQUIRE_UNARY(checkApproxEq(
578 0 : prob.getGradient(), weightFactor * x0 - weightFactor * dc));
579 0 : } else {
580 0 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
581 0 : as<data_t>(0.5) * weightFactor
582 0 : * x0.squaredL2Norm()));
583 :
584 0 : auto gradient = prob.getGradient();
585 0 : DataContainer gradientDirect = weightFactor * x0;
586 0 : for (index_t i = 0; i < gradient.getSize(); ++i)
587 0 : REQUIRE_UNARY(
588 0 : checkApproxEq(gradient[i], gradientDirect[i]));
589 0 : }
590 0 : REQUIRE_EQ(prob.getHessian(), leaf(weightingOp));
591 0 : }
592 2 : } else {
593 2 : if (res.hasOperator()) {
594 0 : if (res.hasDataVector()) {
595 0 : REQUIRE_UNARY(checkApproxEq(
596 0 : prob.evaluate(), as<data_t>(0.5) * scaleFactor * scaleFactor
597 0 : * x0.squaredL2Norm()
598 0 : - scaleFactor * x0.dot(dc)));
599 0 : auto gradient = prob.getGradient();
600 0 : DataContainer gradientDirect =
601 0 : scaleFactor * (scaleFactor * x0) - scaleFactor * dc;
602 0 : for (index_t i = 0; i < gradient.getSize(); ++i)
603 0 : REQUIRE_UNARY(
604 0 : checkApproxEq(gradient[i], gradientDirect[i]));
605 0 : } else {
606 0 : REQUIRE_UNARY(checkApproxEq(
607 0 : prob.evaluate(), as<data_t>(0.5) * scaleFactor * scaleFactor
608 0 : * x0.squaredL2Norm()));
609 0 : auto gradient = prob.getGradient();
610 0 : DataContainer gradientDirect = scaleFactor * (scaleFactor * x0);
611 0 : for (index_t i = 0; i < gradient.getSize(); ++i)
612 0 : REQUIRE_UNARY(
613 0 : checkApproxEq(gradient[i], gradientDirect[i]));
614 0 : }
615 0 : REQUIRE_EQ(prob.getHessian(), leaf(adjoint(scalingOp) * scalingOp));
616 2 : } else {
617 2 : if (res.hasDataVector()) {
618 0 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
619 0 : as<data_t>(0.5) * x0.squaredL2Norm()
620 0 : - x0.dot(dc)));
621 0 : REQUIRE_UNARY(isApprox(prob.getGradient(), x0 - dc));
622 2 : } else {
623 2 : REQUIRE_UNARY(checkApproxEq(
624 2 : prob.evaluate(), as<data_t>(0.5) * x0.squaredL2Norm()));
625 2 : REQUIRE_UNARY(isApprox(prob.getGradient(), x0));
626 2 : }
627 2 : REQUIRE_EQ(prob.getHessian(), leaf(Identity<data_t>{dd}));
628 2 : }
629 2 : }
630 2 : }
631 2 : }
632 2 : }
633 :
634 1664 : for (const auto& regTerm : dataTerms) {
635 1664 : const auto& res =
636 1664 : static_cast<const LinearResidual<data_t>&>(regTerm->getResidual());
637 1664 : const auto isWeighted = is<WeightedL2NormPow2<data_t>>(regTerm.get());
638 1664 : std::string regType =
639 1664 : isWeighted ? "weighted l2-norm regularizer" : "l2-norm regularizer";
640 :
641 1664 : std::string desc =
642 1664 : regType
643 1664 : + (std::string(res.hasOperator() ? " with an operator" : " with no operator")
644 1664 : + ", ")
645 1664 : + (res.hasDataVector() ? "a data vector" : "no data vector");
646 :
647 1664 : WHEN("Converting a Tikhonov problem with no x0 to a quadric problem")
648 1664 : {
649 :
650 2 : INFO("Converting a Tikhonov problem with a ", desc,
651 2 : "and no x0 to a quadric problem");
652 :
653 : // conversion of data terms already tested, fix to 0.5*||Ax-b||^2 for tikhonov
654 : // problem tests
655 2 : auto dataScaleFactor = static_cast<data_t>(5.0);
656 2 : Scaling<data_t> dataScalingOp{dd, dataScaleFactor};
657 :
658 2 : randomData.setRandom();
659 2 : DataContainer<data_t> dataB{dd, randomData};
660 2 : L2NormPow2<data_t> func{LinearResidual<data_t>{dataScalingOp, dataB}};
661 :
662 2 : auto regWeight = static_cast<data_t>(0.01);
663 2 : RegularizationTerm<data_t> reg{regWeight, *regTerm};
664 :
665 2 : Problem<data_t> initialProb{func, reg};
666 :
667 2 : QuadricProblem<data_t> prob{initialProb};
668 :
669 2 : THEN("the clone works correctly")
670 2 : {
671 2 : auto probClone = prob.clone();
672 :
673 2 : REQUIRE_NE(probClone.get(), &prob);
674 2 : REQUIRE_EQ(*probClone, prob);
675 :
676 2 : AND_THEN("the problem behaves as expected")
677 2 : {
678 2 : DataContainer<data_t> zero{dd};
679 2 : zero = 0;
680 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), zero));
681 :
682 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(), 0));
683 :
684 2 : DataContainer<data_t> grad = -dataScaleFactor * dataB;
685 :
686 2 : if (res.hasDataVector()) {
687 0 : DataContainer<data_t> regGrad = static_cast<data_t>(-1.0) * dc;
688 0 : if (isWeighted) {
689 0 : regGrad *= regWeight * weightFactor;
690 0 : }
691 0 : if (res.hasOperator()) {
692 0 : regGrad *= scaleFactor;
693 0 : }
694 0 : if (!isWeighted) {
695 0 : regGrad *= regWeight;
696 0 : }
697 :
698 0 : REQUIRE_UNARY(isApprox(prob.getGradient(), grad + regGrad));
699 2 : } else {
700 2 : REQUIRE_UNARY(isApprox(prob.getGradient(), grad));
701 2 : }
702 :
703 2 : if (isWeighted) {
704 0 : Scaling<data_t> lambdaW{dd, regWeight * weightFactor};
705 :
706 0 : if (res.hasOperator()) {
707 0 : REQUIRE_EQ(prob.getHessian(),
708 0 : leaf(adjoint(dataScalingOp) * dataScalingOp
709 0 : + adjoint(scalingOp) * lambdaW * scalingOp));
710 0 : } else {
711 0 : REQUIRE_EQ(
712 0 : prob.getHessian(),
713 0 : leaf(adjoint(dataScalingOp) * dataScalingOp + lambdaW));
714 0 : }
715 2 : } else {
716 2 : Scaling<data_t> lambdaOp{dd, regWeight};
717 :
718 2 : if (res.hasOperator()) {
719 0 : REQUIRE_EQ(prob.getHessian(),
720 0 : leaf(adjoint(dataScalingOp) * dataScalingOp
721 0 : + lambdaOp * adjoint(scalingOp) * scalingOp));
722 2 : } else {
723 2 : REQUIRE_EQ(
724 2 : prob.getHessian(),
725 2 : leaf(adjoint(dataScalingOp) * dataScalingOp + lambdaOp));
726 2 : }
727 2 : }
728 2 : }
729 2 : }
730 2 : }
731 1664 : WHEN("Converting a Tikhonov problem with x0 to a quadric problem")
732 1664 : {
733 2 : INFO("Converting a Tikhonov problem with a ", desc,
734 2 : "and x0 to a quadric problem");
735 :
736 2 : randomData.setRandom();
737 2 : DataContainer<data_t> x0{dd, randomData};
738 :
739 : // conversion of data terms already tested, fix to 0.5*||Ax-b||^2 for
740 : // tikhonov problem tests
741 2 : auto dataScaleFactor = static_cast<data_t>(5.0);
742 2 : Scaling<data_t> dataScalingOp{dd, dataScaleFactor};
743 :
744 2 : randomData.setRandom();
745 2 : DataContainer<data_t> dataB{dd, randomData};
746 2 : L2NormPow2<data_t> func{LinearResidual<data_t>{dataScalingOp, dataB}};
747 :
748 2 : auto regWeight = static_cast<data_t>(0.01);
749 2 : RegularizationTerm<data_t> reg{regWeight, *regTerm};
750 :
751 2 : Problem<data_t> initialProb{func, reg, x0};
752 :
753 2 : QuadricProblem<data_t> prob{initialProb};
754 :
755 2 : THEN("the clone works correctly")
756 2 : {
757 2 : auto probClone = prob.clone();
758 :
759 2 : REQUIRE_NE(probClone.get(), &prob);
760 2 : REQUIRE_EQ(*probClone, prob);
761 :
762 2 : AND_THEN("the problem behaves as expected")
763 2 : {
764 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
765 :
766 2 : DataContainer<data_t> Ax = dataScaleFactor * (dataScaleFactor * x0);
767 2 : DataContainer<data_t> b = dataScaleFactor * dataB;
768 :
769 2 : if (isWeighted) {
770 0 : Scaling<data_t> lambdaW{dd, regWeight * weightFactor};
771 :
772 0 : if (res.hasOperator()) {
773 0 : Ax += scaleFactor
774 0 : * (regWeight * weightFactor * (scaleFactor * x0));
775 0 : if (res.hasDataVector()) {
776 0 : b += scaleFactor * (regWeight * weightFactor * dc);
777 0 : }
778 0 : REQUIRE_EQ(prob.getHessian(),
779 0 : leaf(adjoint(dataScalingOp) * dataScalingOp
780 0 : + adjoint(scalingOp) * lambdaW * scalingOp));
781 0 : } else {
782 0 : Ax += weightFactor * regWeight * x0;
783 0 : if (res.hasDataVector()) {
784 0 : b += weightFactor * regWeight * dc;
785 0 : }
786 0 : REQUIRE_EQ(
787 0 : prob.getHessian(),
788 0 : leaf(adjoint(dataScalingOp) * dataScalingOp + lambdaW));
789 0 : }
790 2 : } else {
791 2 : Scaling<data_t> lambdaOp{dd, regWeight};
792 :
793 2 : if (res.hasOperator()) {
794 0 : Ax += regWeight * (scaleFactor * (scaleFactor * x0));
795 0 : if (res.hasDataVector()) {
796 0 : b += regWeight * (scaleFactor * dc);
797 0 : }
798 0 : REQUIRE_EQ(prob.getHessian(),
799 0 : leaf(adjoint(dataScalingOp) * dataScalingOp
800 0 : + lambdaOp * adjoint(scalingOp) * scalingOp));
801 2 : } else {
802 2 : Ax += regWeight * x0;
803 2 : if (res.hasDataVector()) {
804 0 : b += regWeight * dc;
805 0 : }
806 2 : REQUIRE_EQ(
807 2 : prob.getHessian(),
808 2 : leaf(adjoint(dataScalingOp) * dataScalingOp + lambdaOp));
809 2 : }
810 2 : }
811 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
812 2 : as<data_t>(0.5) * x0.dot(Ax) - x0.dot(b)));
813 :
814 2 : auto gradient = prob.getGradient();
815 2 : DataContainer gradientDirect = Ax - b;
816 2004 : for (index_t i = 0; i < gradient.getSize(); ++i)
817 2 : REQUIRE_UNARY(checkApproxEq(gradient[i], gradientDirect[i]));
818 2 : }
819 2 : }
820 2 : }
821 1664 : }
822 208 : }
823 26 : }
824 26 : }
825 :
826 : TEST_SUITE_END();
|
