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 32 : TEST_CASE_TEMPLATE("QuadricProblem: Construction with a Quadric functional", TestType, float,
33 : double)
34 : {
35 : using data_t = TestType;
36 :
37 : // eliminate the timing info from console for the tests
38 8 : Logger::setLevel(Logger::LogLevel::WARN);
39 :
40 16 : GIVEN("a Quadric functional")
41 : {
42 16 : IndexVector_t numCoeff(3);
43 8 : numCoeff << 13, 11, 7;
44 16 : VolumeDescriptor dd{numCoeff};
45 :
46 8 : auto scaleFactor = static_cast<data_t>(3.0);
47 16 : Scaling<data_t> scalingOp{dd, scaleFactor};
48 :
49 16 : Eigen::Matrix<data_t, -1, 1> randomData{dd.getNumberOfCoefficients()};
50 8 : randomData.setRandom();
51 16 : DataContainer<data_t> dc{dd, randomData};
52 :
53 16 : Quadric<data_t> quad{scalingOp, dc};
54 :
55 12 : WHEN("constructing a QuadricProblem without x0 from it")
56 : {
57 8 : QuadricProblem<data_t> prob{quad};
58 :
59 6 : THEN("the clone works correctly")
60 : {
61 4 : auto probClone = prob.clone();
62 :
63 2 : REQUIRE_NE(probClone.get(), &prob);
64 2 : REQUIRE_EQ(*probClone, prob);
65 : }
66 :
67 6 : THEN("the problem behaves as expected")
68 : {
69 4 : 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 : }
79 : }
80 :
81 12 : WHEN("constructing a QuadricProblem with x0 from it")
82 : {
83 8 : Eigen::Matrix<data_t, -1, 1> x0Vec{dd.getNumberOfCoefficients()};
84 4 : x0Vec.setRandom();
85 8 : DataContainer<data_t> x0{dd, x0Vec};
86 :
87 8 : QuadricProblem<data_t> prob{quad, x0};
88 :
89 6 : THEN("the clone works correctly")
90 : {
91 4 : auto probClone = prob.clone();
92 :
93 2 : REQUIRE_NE(probClone.get(), &prob);
94 2 : REQUIRE_EQ(*probClone, prob);
95 : }
96 :
97 6 : THEN("the problem behaves as expected")
98 : {
99 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
100 :
101 2 : REQUIRE_UNARY(checkApproxEq(
102 : prob.evaluate(),
103 : static_cast<data_t>(0.5 * scaleFactor) * x0.squaredL2Norm() - x0.dot(dc)));
104 :
105 4 : auto gradient = prob.getGradient();
106 4 : DataContainer gradientDirect = scaleFactor * x0 - dc;
107 2004 : for (index_t i = 0; i < gradient.getSize(); ++i)
108 2002 : REQUIRE_UNARY(checkApproxEq(gradient[i], gradientDirect[i]));
109 :
110 2 : auto hessian = prob.getHessian();
111 2 : REQUIRE_EQ(hessian, leaf(scalingOp));
112 : }
113 : }
114 : }
115 8 : }
116 :
117 32 : TEST_CASE_TEMPLATE("QuadricProblem: with a Quadric functional with spd operator", TestType, float,
118 : double)
119 : {
120 : using data_t = TestType;
121 :
122 : // eliminate the timing info from console for the tests
123 8 : Logger::setLevel(Logger::LogLevel::WARN);
124 :
125 16 : GIVEN("a spd operator")
126 : {
127 16 : IndexVector_t numCoeff(3);
128 8 : numCoeff << 13, 11, 7;
129 16 : VolumeDescriptor dd{numCoeff};
130 :
131 8 : data_t scaleFactor = 3.0;
132 16 : Scaling<data_t> scalingOp{dd, scaleFactor};
133 :
134 16 : Eigen::Matrix<data_t, -1, 1> randomData{dd.getNumberOfCoefficients()};
135 8 : randomData.setRandom();
136 16 : DataContainer<data_t> dc{dd, randomData};
137 :
138 12 : WHEN("constructing a QuadricProblem without x0 from it")
139 : {
140 8 : QuadricProblem<data_t> prob{scalingOp, dc, true};
141 :
142 6 : THEN("the clone works correctly")
143 : {
144 4 : auto probClone = prob.clone();
145 :
146 2 : REQUIRE_NE(probClone.get(), &prob);
147 2 : REQUIRE_EQ(*probClone, prob);
148 : }
149 :
150 6 : THEN("the problem behaves as expected")
151 : {
152 4 : 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 : }
162 : }
163 :
164 12 : WHEN("constructing a QuadricProblem with x0 from it")
165 : {
166 8 : Eigen::Matrix<data_t, -1, 1> x0Vec{dd.getNumberOfCoefficients()};
167 4 : x0Vec.setRandom();
168 8 : DataContainer<data_t> x0{dd, x0Vec};
169 :
170 8 : QuadricProblem<data_t> prob{scalingOp, dc, x0, true};
171 :
172 6 : THEN("the clone works correctly")
173 : {
174 4 : auto probClone = prob.clone();
175 :
176 2 : REQUIRE_NE(probClone.get(), &prob);
177 2 : REQUIRE_EQ(*probClone, prob);
178 : }
179 :
180 6 : THEN("the problem behaves as expected")
181 : {
182 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
183 :
184 2 : REQUIRE_UNARY(
185 : checkApproxEq(prob.evaluate(),
186 : as<data_t>(0.5 * scaleFactor) * x0.squaredL2Norm() - x0.dot(dc)));
187 :
188 4 : auto gradient = prob.getGradient();
189 4 : DataContainer gradientDirect = scaleFactor * x0 - dc;
190 2004 : for (index_t i = 0; i < gradient.getSize(); ++i)
191 2002 : REQUIRE_UNARY(checkApproxEq(gradient[i], gradientDirect[i]));
192 :
193 2 : auto hessian = prob.getHessian();
194 2 : REQUIRE_EQ(hessian, leaf(scalingOp));
195 : }
196 : }
197 : }
198 8 : }
199 :
200 32 : TEST_CASE_TEMPLATE("QuadricProblem: with a Quadric functional with non-spd operator", TestType,
201 : float, double)
202 : {
203 : using data_t = TestType;
204 :
205 : // eliminate the timing info from console for the tests
206 8 : Logger::setLevel(Logger::LogLevel::WARN);
207 :
208 16 : GIVEN("a non-spd operator")
209 : {
210 16 : IndexVector_t numCoeff(3);
211 8 : numCoeff << 13, 11, 7;
212 16 : VolumeDescriptor dd{numCoeff};
213 :
214 8 : data_t scaleFactor = -3.0;
215 16 : Scaling<data_t> scalingOp{dd, scaleFactor};
216 :
217 16 : Eigen::Matrix<data_t, -1, 1> randomData{dd.getNumberOfCoefficients()};
218 8 : randomData.setRandom();
219 16 : DataContainer<data_t> dc{dd, randomData};
220 :
221 12 : WHEN("Constructing a QuadricProblem without x0 from it")
222 : {
223 8 : QuadricProblem<data_t> prob{scalingOp, dc, false};
224 :
225 6 : THEN("the clone works correctly")
226 : {
227 4 : auto probClone = prob.clone();
228 :
229 2 : REQUIRE_NE(probClone.get(), &prob);
230 2 : REQUIRE_EQ(*probClone, prob);
231 : }
232 :
233 6 : THEN("the problem behaves as expected")
234 : {
235 4 : 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 : }
245 : }
246 :
247 12 : WHEN("Constructing a QuadricProblem with x0 from it")
248 : {
249 8 : Eigen::Matrix<data_t, -1, 1> x0Vec{dd.getNumberOfCoefficients()};
250 4 : x0Vec.setRandom();
251 8 : DataContainer<data_t> x0{dd, x0Vec};
252 :
253 8 : QuadricProblem<data_t> prob{scalingOp, dc, x0, false};
254 :
255 6 : THEN("the clone works correctly")
256 : {
257 4 : auto probClone = prob.clone();
258 :
259 2 : REQUIRE_NE(probClone.get(), &prob);
260 2 : REQUIRE_EQ(*probClone, prob);
261 : }
262 :
263 6 : THEN("the problem behaves as expected")
264 : {
265 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
266 :
267 2 : REQUIRE_UNARY(
268 : checkApproxEq(prob.evaluate(),
269 : as<data_t>(0.5 * scaleFactor * scaleFactor) * x0.squaredL2Norm()
270 : - scaleFactor * x0.dot(dc)));
271 :
272 4 : auto gradient = prob.getGradient();
273 4 : DataContainer gradientDirect = scaleFactor * (scaleFactor * x0) - scaleFactor * dc;
274 2004 : for (index_t i = 0; i < gradient.getSize(); ++i)
275 2002 : REQUIRE_UNARY(checkApproxEq(gradient[i], gradientDirect[i]));
276 :
277 2 : auto hessian = prob.getHessian();
278 2 : REQUIRE_EQ(hessian, leaf(adjoint(scalingOp) * scalingOp));
279 : }
280 : }
281 : }
282 8 : }
283 :
284 50 : TEST_CASE_TEMPLATE("QuadricProblem: with a different optimization problems", TestType, float,
285 : double)
286 : {
287 : using data_t = TestType;
288 :
289 : // eliminate the timing info from console for the tests
290 26 : Logger::setLevel(Logger::LogLevel::WARN);
291 :
292 52 : GIVEN("an optimization problem")
293 : {
294 52 : IndexVector_t numCoeff(3);
295 26 : numCoeff << 13, 11, 7;
296 52 : VolumeDescriptor dd{numCoeff};
297 :
298 26 : data_t scaleFactor = -3.0;
299 52 : Scaling<data_t> scalingOp{dd, scaleFactor};
300 :
301 52 : Eigen::Matrix<data_t, -1, 1> randomData{dd.getNumberOfCoefficients()};
302 26 : randomData.setRandom();
303 52 : DataContainer<data_t> dc{dd, randomData};
304 :
305 26 : data_t weightFactor = 2.0;
306 52 : Scaling<data_t> weightingOp{dd, weightFactor};
307 :
308 28 : WHEN("trying to convert a problem with a non-quadric non-wls data term")
309 : {
310 4 : Problem prob{Huber<data_t>{dd}};
311 6 : THEN("an exception is thrown") { REQUIRE_THROWS(QuadricProblem{prob}); }
312 : }
313 :
314 28 : WHEN("trying to convert a problem with a non-quadric non-wls regularization term")
315 : {
316 4 : Problem prob{L2NormPow2<data_t>{dd},
317 : RegularizationTerm{static_cast<data_t>(0.5), Huber<data_t>{dd}}};
318 6 : THEN("an exception is thrown") { REQUIRE_THROWS(QuadricProblem{prob}); }
319 : }
320 :
321 30 : WHEN("converting an optimization problem that has only a quadric data term")
322 : {
323 4 : randomData.setRandom();
324 8 : DataContainer<data_t> x0{dd, randomData};
325 :
326 8 : Problem<data_t> initialProb{Quadric<data_t>{weightingOp, dc}, x0};
327 :
328 8 : QuadricProblem<data_t> prob{initialProb};
329 6 : THEN("the clone works correctly")
330 : {
331 4 : auto probClone = prob.clone();
332 :
333 2 : REQUIRE_NE(probClone.get(), &prob);
334 2 : REQUIRE_EQ(*probClone, prob);
335 : }
336 :
337 6 : THEN("the problem behaves as expected")
338 : {
339 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
340 :
341 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
342 : as<data_t>(0.5 * weightFactor) * x0.squaredL2Norm()
343 : - 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 : }
349 : }
350 :
351 32 : WHEN("calling the conversion constructor on a QadricProblem")
352 : {
353 6 : randomData.setRandom();
354 12 : DataContainer<data_t> x0{dd, randomData};
355 :
356 12 : QuadricProblem<data_t> initialProb{Quadric<data_t>{weightingOp, dc}, x0};
357 :
358 12 : QuadricProblem<data_t> prob{initialProb};
359 8 : THEN("conversion yields the same result as cloning")
360 : {
361 2 : REQUIRE_EQ(*initialProb.clone(), prob);
362 : }
363 :
364 8 : THEN("the clone works correctly")
365 : {
366 4 : auto probClone = prob.clone();
367 :
368 2 : REQUIRE_NE(probClone.get(), &prob);
369 2 : REQUIRE_EQ(*probClone, prob);
370 : }
371 :
372 8 : THEN("the problem behaves as expected")
373 : {
374 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
375 :
376 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
377 : as<data_t>(0.5 * weightFactor) * x0.squaredL2Norm()
378 : - 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 : }
384 : }
385 :
386 26 : randomData.setRandom();
387 52 : DataContainer<data_t> x0{dd, randomData};
388 :
389 52 : Scaling A{dd, static_cast<data_t>(2.0)};
390 26 : randomData.setRandom();
391 52 : DataContainer<data_t> b{dd, randomData};
392 :
393 52 : 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 52 : Scaling L{dd, static_cast<data_t>(0.5)};
400 26 : randomData.setRandom();
401 52 : DataContainer c{dd, randomData};
402 :
403 52 : 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 : "has no operator, but has a vector",
411 : "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 418 : WHEN("Calling the conversion constructor with an OptimizationProblem with only "
416 : "quadric terms")
417 : {
418 4 : INFO("The data term: ", descriptions[i]);
419 4 : INFO("The regularization term: ", descriptions[j]);
420 :
421 4 : RegularizationTerm reg{static_cast<data_t>(0.25), *regFunc[j]};
422 4 : Problem prob{*dataTerms[i], reg, x0};
423 :
424 4 : QuadricProblem converted{prob};
425 :
426 4 : THEN("the problem can be converted and all operations yield the same result as "
427 : "for the initial problem")
428 : {
429 :
430 2 : REQUIRE_UNARY(checkApproxEq(converted.evaluate(), prob.evaluate()));
431 :
432 2 : DataContainer<data_t> gradDiff =
433 : prob.getGradient() - converted.getGradient();
434 2 : REQUIRE_UNARY(checkApproxEq(gradDiff.squaredL2Norm(), 0));
435 :
436 2 : REQUIRE_UNARY(isApprox(prob.getHessian().apply(x0),
437 : converted.getHessian().apply(x0)));
438 : }
439 : }
440 : }
441 : }
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 52 : std::move(std::make_unique<L2NormPow2<data_t>>(LinearResidual<data_t>{scalingOp})));
447 26 : dataTerms.push_back(
448 52 : std::move(std::make_unique<L2NormPow2<data_t>>(LinearResidual<data_t>{dc})));
449 26 : dataTerms.push_back(
450 52 : 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 : LinearResidual<data_t>{scalingOp}, weightingOp)));
454 26 : dataTerms.push_back(std::move(
455 : 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 : LinearResidual<data_t>{scalingOp, dc}, weightingOp)));
458 :
459 234 : 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 416 : std::string probType = isWeighted ? "wls problem" : "ls problem";
463 :
464 832 : std::string desc =
465 : probType
466 208 : + (std::string(res.hasOperator() ? " with an operator" : " with no operator")
467 : + ", ")
468 208 : + (res.hasDataVector() ? "a data vector" : "no data vector");
469 210 : WHEN("Converting a quadric problem with no x0")
470 : {
471 4 : 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 4 : Problem<data_t> initialProb{*dataTerm};
475 :
476 4 : QuadricProblem<data_t> prob{initialProb};
477 :
478 4 : THEN("the clone works correctly")
479 : {
480 4 : auto probClone = prob.clone();
481 :
482 2 : REQUIRE_NE(probClone.get(), &prob);
483 2 : REQUIRE_EQ(*probClone, prob);
484 :
485 4 : AND_THEN("the problem behaves as expected")
486 : {
487 4 : 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 6 : 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 : leaf(adjoint(scalingOp) * weightingOp * scalingOp));
505 : } else {
506 0 : REQUIRE_EQ(prob.getHessian(), leaf(weightingOp));
507 : }
508 : } else {
509 2 : if (res.hasOperator()) {
510 0 : REQUIRE_EQ(prob.getHessian(), leaf(adjoint(scalingOp) * scalingOp));
511 : } else {
512 2 : REQUIRE_EQ(prob.getHessian(), leaf(Identity<data_t>{dd}));
513 : }
514 : }
515 : }
516 : }
517 : }
518 :
519 210 : WHEN("Converting a quadric problem with x0")
520 : {
521 4 : INFO("Converting a ", desc, "and x0 to a quadric problem");
522 2 : randomData.setRandom();
523 4 : DataContainer<data_t> x0{dd, randomData};
524 :
525 : // use OptimizationProblem istead of WLSProblem as it allows for more general
526 : // formulations
527 4 : Problem<data_t> initialProb{*dataTerm, x0};
528 :
529 4 : QuadricProblem<data_t> prob{initialProb};
530 :
531 4 : THEN("the clone works correctly")
532 : {
533 4 : auto probClone = prob.clone();
534 :
535 2 : REQUIRE_NE(probClone.get(), &prob);
536 2 : REQUIRE_EQ(*probClone, prob);
537 :
538 4 : AND_THEN("the problem behaves as expected")
539 : {
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 : checkApproxEq(prob.evaluate(),
547 : (as<data_t>(0.5) * scaleFactor * scaleFactor
548 : * weightFactor * x0.squaredL2Norm()
549 : - scaleFactor * weightFactor * x0.dot(dc))));
550 0 : auto gradient = prob.getGradient();
551 0 : DataContainer gradientDirect =
552 : scaleFactor * (weightFactor * (scaleFactor * x0))
553 : - scaleFactor * (weightFactor * dc);
554 0 : for (index_t i = 0; i < gradient.getSize(); ++i)
555 0 : REQUIRE_UNARY(
556 : checkApproxEq(gradient[i], gradientDirect[i]));
557 : } else {
558 0 : REQUIRE_UNARY(checkApproxEq(
559 : prob.evaluate(), as<data_t>(0.5) * scaleFactor * scaleFactor
560 : * weightFactor * x0.squaredL2Norm()));
561 :
562 0 : auto gradient = prob.getGradient();
563 0 : DataContainer gradientDirect =
564 : scaleFactor * (weightFactor * (scaleFactor * x0));
565 0 : for (index_t i = 0; i < gradient.getSize(); ++i)
566 0 : REQUIRE_UNARY(
567 : checkApproxEq(gradient[i], gradientDirect[i]));
568 : }
569 0 : REQUIRE_EQ(prob.getHessian(),
570 : leaf(adjoint(scalingOp) * weightingOp * scalingOp));
571 : } else {
572 0 : if (res.hasDataVector()) {
573 0 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
574 : as<data_t>(0.5) * weightFactor
575 : * x0.squaredL2Norm()
576 : - weightFactor * x0.dot(dc)));
577 0 : REQUIRE_UNARY(checkApproxEq(
578 : prob.getGradient(), weightFactor * x0 - weightFactor * dc));
579 : } else {
580 0 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
581 : as<data_t>(0.5) * weightFactor
582 : * 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 : checkApproxEq(gradient[i], gradientDirect[i]));
589 : }
590 0 : REQUIRE_EQ(prob.getHessian(), leaf(weightingOp));
591 : }
592 : } else {
593 2 : if (res.hasOperator()) {
594 0 : if (res.hasDataVector()) {
595 0 : REQUIRE_UNARY(checkApproxEq(
596 : prob.evaluate(), as<data_t>(0.5) * scaleFactor * scaleFactor
597 : * x0.squaredL2Norm()
598 : - scaleFactor * x0.dot(dc)));
599 0 : auto gradient = prob.getGradient();
600 0 : DataContainer gradientDirect =
601 : scaleFactor * (scaleFactor * x0) - scaleFactor * dc;
602 0 : for (index_t i = 0; i < gradient.getSize(); ++i)
603 0 : REQUIRE_UNARY(
604 : checkApproxEq(gradient[i], gradientDirect[i]));
605 : } else {
606 0 : REQUIRE_UNARY(checkApproxEq(
607 : prob.evaluate(), as<data_t>(0.5) * scaleFactor * scaleFactor
608 : * 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 : checkApproxEq(gradient[i], gradientDirect[i]));
614 : }
615 0 : REQUIRE_EQ(prob.getHessian(), leaf(adjoint(scalingOp) * scalingOp));
616 : } else {
617 2 : if (res.hasDataVector()) {
618 0 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
619 : as<data_t>(0.5) * x0.squaredL2Norm()
620 : - x0.dot(dc)));
621 0 : REQUIRE_UNARY(isApprox(prob.getGradient(), x0 - dc));
622 : } else {
623 2 : REQUIRE_UNARY(checkApproxEq(
624 : prob.evaluate(), as<data_t>(0.5) * x0.squaredL2Norm()));
625 2 : REQUIRE_UNARY(isApprox(prob.getGradient(), x0));
626 : }
627 2 : REQUIRE_EQ(prob.getHessian(), leaf(Identity<data_t>{dd}));
628 : }
629 : }
630 : }
631 : }
632 : }
633 :
634 1872 : for (const auto& regTerm : dataTerms) {
635 : 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 3328 : std::string regType =
639 : isWeighted ? "weighted l2-norm regularizer" : "l2-norm regularizer";
640 :
641 6656 : std::string desc =
642 : regType
643 1664 : + (std::string(res.hasOperator() ? " with an operator" : " with no operator")
644 : + ", ")
645 1664 : + (res.hasDataVector() ? "a data vector" : "no data vector");
646 :
647 1666 : WHEN("Converting a Tikhonov problem with no x0 to a quadric problem")
648 : {
649 :
650 4 : INFO("Converting a Tikhonov problem with a ", desc,
651 : "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 4 : Scaling<data_t> dataScalingOp{dd, dataScaleFactor};
657 :
658 2 : randomData.setRandom();
659 4 : DataContainer<data_t> dataB{dd, randomData};
660 4 : L2NormPow2<data_t> func{LinearResidual<data_t>{dataScalingOp, dataB}};
661 :
662 2 : auto regWeight = static_cast<data_t>(0.01);
663 4 : RegularizationTerm<data_t> reg{regWeight, *regTerm};
664 :
665 4 : Problem<data_t> initialProb{func, reg};
666 :
667 4 : QuadricProblem<data_t> prob{initialProb};
668 :
669 4 : THEN("the clone works correctly")
670 : {
671 4 : auto probClone = prob.clone();
672 :
673 2 : REQUIRE_NE(probClone.get(), &prob);
674 2 : REQUIRE_EQ(*probClone, prob);
675 :
676 4 : AND_THEN("the problem behaves as expected")
677 : {
678 4 : 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 4 : 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 : }
691 0 : if (res.hasOperator()) {
692 0 : regGrad *= scaleFactor;
693 : }
694 0 : if (!isWeighted) {
695 0 : regGrad *= regWeight;
696 : }
697 :
698 0 : REQUIRE_UNARY(isApprox(prob.getGradient(), grad + regGrad));
699 : } else {
700 2 : REQUIRE_UNARY(isApprox(prob.getGradient(), grad));
701 : }
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 : leaf(adjoint(dataScalingOp) * dataScalingOp
709 : + adjoint(scalingOp) * lambdaW * scalingOp));
710 : } else {
711 0 : REQUIRE_EQ(
712 : prob.getHessian(),
713 : leaf(adjoint(dataScalingOp) * dataScalingOp + lambdaW));
714 : }
715 : } else {
716 4 : Scaling<data_t> lambdaOp{dd, regWeight};
717 :
718 2 : if (res.hasOperator()) {
719 0 : REQUIRE_EQ(prob.getHessian(),
720 : leaf(adjoint(dataScalingOp) * dataScalingOp
721 : + lambdaOp * adjoint(scalingOp) * scalingOp));
722 : } else {
723 2 : REQUIRE_EQ(
724 : prob.getHessian(),
725 : leaf(adjoint(dataScalingOp) * dataScalingOp + lambdaOp));
726 : }
727 : }
728 : }
729 : }
730 : }
731 1666 : WHEN("Converting a Tikhonov problem with x0 to a quadric problem")
732 : {
733 4 : INFO("Converting a Tikhonov problem with a ", desc,
734 : "and x0 to a quadric problem");
735 :
736 2 : randomData.setRandom();
737 4 : 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 4 : Scaling<data_t> dataScalingOp{dd, dataScaleFactor};
743 :
744 2 : randomData.setRandom();
745 4 : DataContainer<data_t> dataB{dd, randomData};
746 4 : L2NormPow2<data_t> func{LinearResidual<data_t>{dataScalingOp, dataB}};
747 :
748 2 : auto regWeight = static_cast<data_t>(0.01);
749 4 : RegularizationTerm<data_t> reg{regWeight, *regTerm};
750 :
751 4 : Problem<data_t> initialProb{func, reg, x0};
752 :
753 4 : QuadricProblem<data_t> prob{initialProb};
754 :
755 4 : THEN("the clone works correctly")
756 : {
757 4 : auto probClone = prob.clone();
758 :
759 2 : REQUIRE_NE(probClone.get(), &prob);
760 2 : REQUIRE_EQ(*probClone, prob);
761 :
762 4 : AND_THEN("the problem behaves as expected")
763 : {
764 2 : REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), x0));
765 :
766 4 : DataContainer<data_t> Ax = dataScaleFactor * (dataScaleFactor * x0);
767 4 : 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 : }
778 0 : REQUIRE_EQ(prob.getHessian(),
779 : leaf(adjoint(dataScalingOp) * dataScalingOp
780 : + adjoint(scalingOp) * lambdaW * scalingOp));
781 : } else {
782 0 : Ax += weightFactor * regWeight * x0;
783 0 : if (res.hasDataVector()) {
784 0 : b += weightFactor * regWeight * dc;
785 : }
786 0 : REQUIRE_EQ(
787 : prob.getHessian(),
788 : leaf(adjoint(dataScalingOp) * dataScalingOp + lambdaW));
789 : }
790 : } else {
791 4 : 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 : }
798 0 : REQUIRE_EQ(prob.getHessian(),
799 : leaf(adjoint(dataScalingOp) * dataScalingOp
800 : + lambdaOp * adjoint(scalingOp) * scalingOp));
801 : } else {
802 2 : Ax += regWeight * x0;
803 2 : if (res.hasDataVector()) {
804 0 : b += regWeight * dc;
805 : }
806 2 : REQUIRE_EQ(
807 : prob.getHessian(),
808 : leaf(adjoint(dataScalingOp) * dataScalingOp + lambdaOp));
809 : }
810 : }
811 2 : REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
812 : as<data_t>(0.5) * x0.dot(Ax) - x0.dot(b)));
813 :
814 4 : auto gradient = prob.getGradient();
815 4 : DataContainer gradientDirect = Ax - b;
816 2004 : for (index_t i = 0; i < gradient.getSize(); ++i)
817 2002 : REQUIRE_UNARY(checkApproxEq(gradient[i], gradientDirect[i]));
818 : }
819 : }
820 : }
821 : }
822 : }
823 : }
824 26 : }
825 :
826 : TEST_SUITE_END();
|
