Line data Source code
1 : /*
2 : * File: TestQuadratic.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on Jul 9, 2015, 4:14:39 AM
6 : *
7 : * ForBES is free software: you can redistribute it and/or modify
8 : * it under the terms of the GNU Lesser General Public License as published by
9 : * the Free Software Foundation, either version 3 of the License, or
10 : * (at your option) any later version.
11 : *
12 : * ForBES is distributed in the hope that it will be useful,
13 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : * GNU Lesser General Public License for more details.
16 : *
17 : * You should have received a copy of the GNU Lesser General Public License
18 : * along with ForBES. If not, see <http://www.gnu.org/licenses/>.
19 : */
20 :
21 : #include <complex>
22 :
23 : #include "TestQuadratic.h"
24 :
25 :
26 : const static double MAT1[16] = {
27 : 7, 2, -2, -1,
28 : 2, 3, 0, -1,
29 : -2, 0, 3, -1,
30 : -1, -1, -1, 1
31 : };
32 : const static double MAT2[16] = {
33 : 16.0, 2.0, 3.0, 13.0,
34 : 5.0, 11.0, 10.0, 8.0,
35 : 9.0, 7.0, 6.0, 12.0,
36 : 4.0, 14.0, 15.0, 1.0
37 : };
38 :
39 1 : CPPUNIT_TEST_SUITE_REGISTRATION(TestQuadratic);
40 :
41 18 : TestQuadratic::TestQuadratic() {
42 :
43 18 : }
44 :
45 36 : TestQuadratic::~TestQuadratic() {
46 36 : }
47 :
48 18 : void TestQuadratic::setUp() {
49 18 : }
50 :
51 18 : void TestQuadratic::tearDown() {
52 18 : }
53 :
54 1 : void TestQuadratic::testQuadratic() {
55 : const double * Qdata;
56 1 : Qdata = MAT2;
57 1 : double xdata[4] = {-1.0, 1.0, 1.0, 1.0};
58 :
59 1 : Matrix Q = Matrix(4, 4, Qdata);
60 2 : Matrix x = Matrix(4, 1, xdata);
61 1 : Function *quad = NULL;
62 1 : _ASSERT_OK(quad = new Quadratic(Q));
63 1 : _ASSERT_NEQ(quad, NULL);
64 :
65 : double f;
66 1 : _ASSERT(quad->category().defines_f());
67 1 : int info = quad -> call(x, f);
68 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, info);
69 :
70 1 : if (quad != NULL) {
71 1 : delete quad;
72 1 : }
73 1 : }
74 :
75 1 : void TestQuadratic::testQuadratic2() {
76 : /* Test the empty constructor */
77 1 : size_t n = 8;
78 1 : const double tol = 1e-7;
79 1 : Function *F = new Quadratic(); /* Q = I, q = 0. */
80 :
81 1 : Matrix x = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0, Matrix::MATRIX_DENSE);
82 1 : x.set(0, 0, 666);
83 2 : Matrix grad;
84 :
85 1 : double fval = 0.0;
86 1 : _ASSERT(F->category().defines_f());
87 1 : _ASSERT(F->category().defines_grad());
88 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, F->call(x, fval, grad));
89 :
90 1 : const double fval_expected = static_cast<Matrix> (x * x).get(0, 0) / 2;
91 1 : _ASSERT_NUM_EQ(fval_expected, fval, tol);
92 :
93 2 : Matrix q(n, 1);
94 9 : for (size_t i = 0; i < n; ++i) {
95 8 : _ASSERT_OK(q.set(i, 0, i + 1));
96 : }
97 1 : _ASSERT_OK(static_cast<Quadratic*> (F)->setq(q));
98 :
99 1 : double fval_prev = fval;
100 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, F->call(x, fval, grad));
101 1 : _ASSERT_NUM_EQ(fval_prev + (q * x).get(0, 0), fval, 1e-6);
102 :
103 2 : _ASSERT_OK(delete F);
104 1 : }
105 :
106 1 : void TestQuadratic::testQuadratic3() {
107 1 : size_t n = 10;
108 1 : Matrix Q(n, n, Matrix::MATRIX_DENSE);
109 10 : for (size_t i = 0; i < n - 1; ++i) {
110 9 : Q.set(i, i, 1.5);
111 : }
112 1 : Quadratic *F = new Quadratic(Q);
113 2 : Matrix x(n, 1);
114 : double fval;
115 1 : _ASSERT(F->category().defines_conjugate());
116 1 : int status = F->callConj(x, fval);
117 1 : _ASSERT_EQ(ForBESUtils::STATUS_NUMERICAL_PROBLEMS, status);
118 :
119 1 : Q.set(n - 1, n - 1, 1.5);
120 1 : Q.set(1, 0, 0.1);
121 1 : Q.set(0, 1, 0.1);
122 :
123 1 : F->setQ(Q);
124 1 : _ASSERT_OK(status = F->callConj(x, fval));
125 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, status);
126 :
127 2 : Matrix q(n, 1);
128 11 : for (size_t i = 0; i < n; ++i) {
129 10 : q.set(i, 0, i + 1);
130 10 : x.set(i, 0, 2 * i + 1);
131 : }
132 :
133 :
134 2 : Matrix grad;
135 1 : static_cast<Quadratic*> (F)->setq(q);
136 :
137 1 : _ASSERT_OK(status = F->callConj(x, fval, grad));
138 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, status);
139 :
140 1 : const double fval_exp = 190.002976190476;
141 1 : const double tol = 1e-8;
142 :
143 1 : _ASSERT_NUM_EQ(fval_exp, fval, tol);
144 :
145 1 : double grad0_exp = -0.0446428571428572;
146 1 : _ASSERT_NUM_EQ(grad0_exp, grad.get(0, 0), tol);
147 :
148 2 : _ASSERT_OK(delete F);
149 1 : }
150 :
151 1 : void TestQuadratic::testCallProx() {
152 1 : Function *F = new Quadratic(); // F(x) = (1/2)x'x
153 1 : Matrix x(4, 1);
154 5 : for (size_t i = 0; i < 4; i++) {
155 4 : x[i] = 0.34 + std::sqrt(i+0.3);
156 : }
157 :
158 1 : const double gamma = 0.737;
159 1 : const double tol = 1e-7;
160 :
161 2 : Matrix prox;
162 :
163 1 : _ASSERT_NOT(F->category().defines_prox());
164 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, F->callProx(x, gamma, prox));
165 1 : _ASSERT_NUM_EQ(0.511066527061120, prox[0], tol);
166 1 : _ASSERT_NUM_EQ(0.852144746746769, prox[1], tol);
167 1 : _ASSERT_NUM_EQ(1.068840005072142, prox[2], tol);
168 1 : _ASSERT_NUM_EQ(1.241560283510935, prox[3], tol);
169 :
170 2 : delete F;
171 1 : }
172 :
173 1 : void TestQuadratic::testCallProx2() {
174 1 : const size_t n = 3;
175 : double data_Q[9] = {
176 : 0.29507822324818300, 0.00199205562668041, 0.00423436821643449,
177 : 0.00199205562668041, 0.24650478030564038, 0.00233990839749903,
178 : 0.00423436821643449, 0.00233990839749903, 0.27841699644617668
179 1 : };
180 1 : Matrix Q(n, n, data_Q);
181 :
182 1 : double data_q[n * n] = {0.1, 0.2, -0.3};
183 2 : Matrix q(n, 1, data_q);
184 :
185 1 : const double gamma = 0.02;
186 :
187 1 : double data_v[n] = {-1.0, 2.0, -0.5};
188 2 : Matrix v(n, 1, data_v);
189 :
190 : double data_prox_exp[] = {
191 : -0.996158636187904,
192 : 1.986270176873906,
193 : -0.491273016601634
194 1 : };
195 2 : Matrix prox_exp(n, 1, data_prox_exp);
196 :
197 2 : Quadratic F(Q, q);
198 2 : Matrix prox(n, 1);
199 1 : F.callProx(v, gamma, prox);
200 :
201 2 : _ASSERT_EQ(prox_exp, prox);
202 1 : }
203 :
204 1 : void TestQuadratic::testCall2() {
205 1 : Function * quad = new Quadratic();
206 1 : Matrix x(4, 1);
207 5 : for (size_t i = 0; i < 4; i++) {
208 4 : x[i] = 0.1 * (i + 1.0);
209 : }
210 : double f;
211 1 : int status = quad->call(x, f);
212 1 : _ASSERT(ForBESUtils::is_status_ok(status));
213 1 : _ASSERT_NUM_EQ(0.15, f, 1e-6);
214 1 : delete quad;
215 1 : }
216 :
217 1 : void TestQuadratic::testCall() {
218 : const double * Qdata;
219 1 : Qdata = MAT2;
220 1 : double qdata[4] = {2.0, 3.0, 4.0, 5.0};
221 1 : double xdata[4] = {-1.0, 1.0, 1.0, 1.0};
222 :
223 1 : Matrix Q = Matrix(4, 4, Qdata);
224 2 : Matrix q = Matrix(4, 1, qdata);
225 2 : Matrix x = Matrix(4, 1, xdata);
226 :
227 2 : Quadratic quadratic(Q, q);
228 1 : double f = -999.0;
229 1 : int status = quadratic.call(x, f);
230 :
231 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, status);
232 1 : _ASSERT_EQ(42.0, f);
233 :
234 : /* Second part */
235 2 : Quadratic quadratic2(Q);
236 1 : status = quadratic2.call(x, f);
237 1 : _ASSERT_EQ(0, status);
238 2 : _ASSERT_EQ(32.0, f);
239 1 : }
240 :
241 1 : void TestQuadratic::testCallWithGradient() {
242 1 : const size_t n = 4;
243 : const double * Qdata;
244 1 : Qdata = MAT1;
245 1 : double xdata[4] = {-1.0, 1.0, 1.0, 1.0};
246 1 : const double expected[4] = {-8, 0, 4, 0};
247 :
248 1 : Matrix Q = Matrix(n, n, Qdata);
249 2 : Matrix x = Matrix(n, 1, xdata);
250 :
251 1 : Function * quad = new Quadratic(Q);
252 1 : double f = -999.0;
253 2 : Matrix grad;
254 1 : _ASSERT(quad->category().defines_f());
255 1 : _ASSERT(quad->category().defines_grad());
256 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, quad -> call(x, f, grad));
257 :
258 5 : for (size_t i = 0; i < n; i++) {
259 4 : _ASSERT_EQ(expected[i], grad[i]);
260 : }
261 :
262 2 : delete quad;
263 :
264 1 : }
265 :
266 1 : void TestQuadratic::testCallConj() {
267 1 : const size_t n = 4;
268 : const double * Qdata;
269 1 : Qdata = MAT1;
270 :
271 1 : double qdata[4] = {2.0, 3.0, 4.0, 5.0};
272 1 : double xdata[4] = {-1.0, 1.0, 1.0, 1.0};
273 :
274 1 : Matrix Q = Matrix(n, n, Qdata);
275 2 : Matrix q = Matrix(n, 1, qdata);
276 2 : Matrix x = Matrix(n, 1, xdata);
277 :
278 2 : Quadratic quadratic(Q, q);
279 : double fstar;
280 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, quadratic.callConj(x, fstar));
281 :
282 :
283 1 : double expected = 421.0;
284 1 : const double rel_tol = 1e-5;
285 1 : _ASSERT(std::fabs(expected - fstar) / expected < rel_tol);
286 :
287 5 : for (size_t i = 0; i < n; i++) {
288 4 : x[i] = 10 + 2 * i;
289 : }
290 :
291 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, quadratic.callConj(x, fstar));
292 1 : expected = 3722;
293 1 : _ASSERT(std::fabs(expected - fstar) / expected < rel_tol);
294 :
295 2 : Matrix grad;
296 :
297 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, quadratic.callConj(x, fstar, grad));
298 :
299 2 : Matrix grad_expected(n, 1);
300 1 : grad_expected.set(0, 0, 45.5);
301 1 : grad_expected.set(1, 0, 35.5);
302 1 : grad_expected.set(2, 0, 96.5);
303 1 : grad_expected.set(3, 0, 188.5);
304 :
305 1 : const double tol = 1e-5;
306 5 : for (size_t i = 0; i < n; i++) {
307 4 : _ASSERT_NUM_EQ(grad_expected.get(i, 0), grad.get(i, 0), tol);
308 1 : }
309 :
310 1 : }
311 :
312 1 : void TestQuadratic::testCallConj2() {
313 1 : const size_t n = 4;
314 :
315 1 : double qdata[4] = {2.0, 3.0, 4.0, 5.0};
316 1 : double xdata[4] = {-1.0, 1.0, 1.0, 1.0};
317 :
318 1 : Quadratic *F = new Quadratic();
319 1 : Matrix q(n, 1, qdata);
320 2 : Matrix x(n, 1, xdata);
321 :
322 1 : _ASSERT_OK(F->setq(q));
323 : double fstar;
324 1 : double fstar_exp = 38.0;
325 1 : const double tol = 1e-6;
326 :
327 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, F->callConj(x, fstar));
328 1 : _ASSERT_NUM_EQ(fstar_exp, fstar, tol);
329 1 : _ASSERT_OK(delete F);
330 :
331 1 : double alpha = 2.3;
332 2 : Matrix Eye = MatrixFactory::MakeIdentity(n, alpha); /* Here Eye is a diagonal matrix */
333 1 : Function *F2 = new Quadratic(Eye, q);
334 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, F2->callConj(x, fstar));
335 1 : _ASSERT_NUM_EQ(fstar_exp / alpha, fstar, tol);
336 2 : _ASSERT_OK(delete F2);
337 1 : }
338 :
339 1 : void TestQuadratic::testCallDiagonalMatrix() {
340 1 : const size_t n = 4;
341 1 : Matrix Q(n, n, Matrix::MATRIX_DIAGONAL);
342 5 : for (size_t i = 0; i < n; i++) {
343 4 : Q[i] = i + 2.0;
344 : }
345 :
346 1 : Function *f = new Quadratic(Q);
347 :
348 2 : Matrix x(n, 1);
349 5 : for (size_t i = 0; i < n; i++) {
350 4 : x[i] = 2 * i + 1.0;
351 : }
352 :
353 : double val;
354 1 : f -> call(x, val);
355 :
356 1 : double tol = 1e-8;
357 1 : _ASSERT_NUM_EQ(187.0, val, tol);
358 2 : _ASSERT_OK(delete f);
359 1 : }
360 :
361 1 : void TestQuadratic::testCallSparse() {
362 : /*
363 : * Q is sparse, x is dense
364 : */
365 1 : size_t n = 10;
366 1 : size_t nnz_Q = 20;
367 1 : Matrix Qsp = MatrixFactory::MakeRandomSparse(n, n, nnz_Q, 0.0, 1.0);
368 :
369 1 : Function *F = new Quadratic(Qsp);
370 :
371 2 : Matrix x = MatrixFactory::MakeRandomMatrix(n, 1, 3.0, 1.5, Matrix::MATRIX_DENSE);
372 1 : double fval = -1;
373 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, F->call(x, fval));
374 1 : _ASSERT(fval > 0);
375 :
376 1 : double f_exp = Qsp.quad(x);
377 1 : const double tol = 1e-8;
378 1 : _ASSERT_NUM_EQ(f_exp, fval, tol);
379 :
380 2 : _ASSERT_OK(delete F);
381 :
382 1 : }
383 :
384 1 : void TestQuadratic::testCallSparse2() {
385 : /*
386 : * BOTH Q and x are sparse
387 : */
388 1 : size_t n = 10;
389 1 : size_t nnz_Q = 20;
390 1 : size_t nnz_x = 6;
391 1 : Matrix Qsp = MatrixFactory::MakeRandomSparse(n, n, nnz_Q, 0.0, 1.0);
392 :
393 1 : Function *F = new Quadratic(Qsp);
394 :
395 2 : Matrix x = MatrixFactory::MakeRandomSparse(n, 1, nnz_x, 0.0, 1.0);
396 1 : double fval = -1;
397 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, F->call(x, fval));
398 1 : _ASSERT(fval > 0);
399 :
400 1 : double f_exp = Qsp.quad(x);
401 1 : const double tol = 1e-8;
402 1 : _ASSERT_NUM_EQ(f_exp, fval, tol);
403 :
404 2 : _ASSERT_OK(delete F);
405 1 : }
406 :
407 1 : void TestQuadratic::testCallConjSparse() {
408 1 : size_t n = 5;
409 1 : size_t nnz_Q = 2 * n - 1;
410 1 : double fstar = -1;
411 1 : const double tol = 1e-8;
412 1 : const double fstar_exp = 5.13142364727941;
413 : int status;
414 : Function *F;
415 :
416 1 : Matrix Qsp = MatrixFactory::MakeSparseSymmetric(n, nnz_Q);
417 2 : Matrix x(n, 1);
418 :
419 6 : for (size_t i = 0; i < n; i++) {
420 5 : Qsp.set(i, i, 10.0);
421 5 : x.set(i, 0, i + 1);
422 : }
423 5 : for (size_t i = 1; i < n; i++) { /* Set the LT part */
424 4 : Qsp.set(i, i - 1, 0.5);
425 : }
426 :
427 1 : F = new Quadratic(Qsp);
428 1 : _ASSERT_OK(status = F->callConj(x, fstar));
429 1 : _ASSERT_NEQ(-1, fstar);
430 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, status);
431 :
432 1 : _ASSERT_NUM_EQ(fstar_exp, fstar, tol);
433 :
434 2 : _ASSERT_OK(delete F);
435 :
436 1 : }
437 :
438 1 : void TestQuadratic::testCallSparse3() {
439 1 : const size_t n = 10;
440 1 : const double tol = 1e-7;
441 1 : Matrix Q = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 2.0, Matrix::MATRIX_SYMMETRIC);
442 2 : Matrix q = MatrixFactory::MakeSparse(n, 1, 0, Matrix::SPARSE_UNSYMMETRIC); /* q = [] */
443 2 : Matrix x = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0, Matrix::MATRIX_DENSE);
444 :
445 1 : Function *F = new Quadratic(Q, q);
446 :
447 : double f;
448 1 : int status = F->call(x, f);
449 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, status);
450 1 : double f_exp = Q.quad(x);
451 1 : _ASSERT_NUM_EQ(f_exp, f, tol);
452 :
453 1 : q = MatrixFactory::MakeSparse(n, 1, 1, Matrix::SPARSE_UNSYMMETRIC);
454 1 : q.set(0, 0, 5.5);
455 :
456 :
457 1 : status = F->call(x, f);
458 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, status);
459 1 : f_exp = Q.quad(x, q);
460 1 : _ASSERT_NUM_EQ(f_exp, f, tol);
461 :
462 2 : _ASSERT_OK(delete F);
463 :
464 1 : }
465 :
466 1 : void TestQuadratic::testHessian() {
467 1 : const size_t n = 4;
468 :
469 1 : Matrix Q = MatrixFactory::MakeRandomMatrix(n, n, 5.0, 1.0);
470 2 : Matrix x = MatrixFactory::MakeRandomMatrix(n, 1, 1.0, 2.0);
471 2 : Matrix d = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0);
472 2 : Matrix Hd = Matrix(n, 1);
473 2 : Matrix Hd_reference = Matrix(n, 1);
474 :
475 1 : Function * quad = new Quadratic(Q);
476 1 : int status = Matrix::mult(Hd_reference, 1.0, Q, d, 0.0);
477 1 : _ASSERT(ForBESUtils::is_status_ok(status));
478 1 : status = quad->hessianProduct(x, d, Hd);
479 1 : _ASSERT(ForBESUtils::is_status_ok(status));
480 :
481 1 : _ASSERT_EQ(Hd, Hd_reference);
482 :
483 2 : delete quad;
484 1 : }
485 :
486 1 : void TestQuadratic::testHessianSparse() {
487 1 : const size_t n = 10;
488 1 : const size_t nnz = 45;
489 :
490 1 : Matrix Q = MatrixFactory::MakeRandomSparse(n, n, nnz, 5.0, 1.0);
491 2 : Matrix x = MatrixFactory::MakeRandomMatrix(n, 1, 1.0, 2.0);
492 2 : Matrix d = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0);
493 2 : Matrix Hd = Matrix(n, 1);
494 2 : Matrix Hd_reference = Matrix(n, 1);
495 :
496 1 : Function * quad = new Quadratic(Q);
497 1 : Hd_reference = Q*d;
498 1 : int status = quad->hessianProduct(x, d, Hd);
499 1 : _ASSERT(ForBESUtils::is_status_ok(status));
500 1 : _ASSERT_EQ(Hd, Hd_reference);
501 :
502 2 : delete quad;
503 1 : }
504 :
505 1 : void TestQuadratic::testApproximateHessian() {
506 :
507 1 : const size_t n = 4;
508 : const double * Qdata;
509 1 : Matrix grad_x;
510 : double f;
511 1 : Qdata = MAT1;
512 :
513 1 : double qdata[4] = {2.0, 3.0, 4.0, 5.0};
514 1 : double xdata[4] = {-1.0, 1.0, 1.0, 1.0};
515 :
516 2 : Matrix Q = Matrix(n, n, Qdata);
517 2 : Matrix q = Matrix(n, 1, qdata);
518 :
519 :
520 2 : Quadratic quadratic(Q, q);
521 :
522 2 : Matrix x = Matrix(n, 1, xdata);
523 2 : Matrix z = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 2.0);
524 1 : const double epsilon = 1e-8;
525 1 : const double one_over_epsilon = 1e8;
526 :
527 : /* t = t + εz */
528 2 : Matrix t(x);
529 1 : Matrix::add(t, epsilon, z, 1.0);
530 :
531 :
532 2 : Matrix Hz(n, 1);
533 1 : quadratic.call(t, f, Hz); /* Hz = nabla f(t) */
534 1 : quadratic.call(x, f, grad_x); /* grad_x = nabla f(x) */
535 1 : Hz -= grad_x; /* Hz = nabla f(t) - nabla f(x) */
536 1 : Hz *= one_over_epsilon;
537 :
538 2 : Matrix Hz_copy = Hz;
539 :
540 1 : quadratic.hessianProduct(x, z, Hz);
541 :
542 5 : for (size_t i = 0; i < n; i++) {
543 4 : _ASSERT_NUM_EQ(Hz_copy[i], Hz[i], 1e-6);
544 1 : }
545 4 : }
546 :
547 :
|
