1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003 Ferdinando Ametrano
5 Copyright (C) 2007, 2008 Klaus Spanderen
6 Copyright (C) 2007 Neil Firth
7 Copyright (C) 2016 Peter Caspers
8
9 This file is part of QuantLib, a free-software/open-source library
10 for financial quantitative analysts and developers - http://quantlib.org/
11
12 QuantLib is free software: you can redistribute it and/or modify it
13 under the terms of the QuantLib license. You should have received a
14 copy of the license along with this program; if not, please email
15 <quantlib-dev@lists.sf.net>. The license is also available online at
16 <http://quantlib.org/license.shtml>.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the license for more details.
21*/
22
23#include "matrices.hpp"
24#include "utilities.hpp"
25#include <ql/experimental/math/moorepenroseinverse.hpp>
26#include <ql/math/matrix.hpp>
27#include <ql/math/matrixutilities/basisincompleteordered.hpp>
28#include <ql/math/matrixutilities/bicgstab.hpp>
29#include <ql/math/matrixutilities/choleskydecomposition.hpp>
30#include <ql/math/matrixutilities/gmres.hpp>
31#include <ql/math/matrixutilities/pseudosqrt.hpp>
32#include <ql/math/matrixutilities/qrdecomposition.hpp>
33#include <ql/math/matrixutilities/svd.hpp>
34#include <ql/math/matrixutilities/symmetricschurdecomposition.hpp>
35#include <ql/math/matrixutilities/sparsematrix.hpp>
36#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
37#include <cmath>
38#include <iterator>
39#include <utility>
40#include <numeric>
41
42using namespace QuantLib;
43using namespace boost::unit_test_framework;
44
45using std::fabs;
46
47#ifdef __cpp_concepts
48static_assert(std::random_access_iterator<Matrix::column_iterator>);
49static_assert(std::random_access_iterator<Matrix::const_column_iterator>);
50static_assert(std::random_access_iterator<Matrix::reverse_column_iterator>);
51static_assert(std::random_access_iterator<Matrix::const_reverse_column_iterator>);
52#endif
53
54namespace matrices_test {
55
56 Size N;
57 Matrix M1, M2, M3, M4, M5, M6, M7, I;
58
59 Real norm(const Array& v) {
60 return std::sqrt(x: DotProduct(v1: v,v2: v));
61 }
62
63 Real norm(const Matrix& m) {
64 Real sum = 0.0;
65 for (Size i=0; i<m.rows(); i++)
66 for (Size j=0; j<m.columns(); j++)
67 sum += m[i][j]*m[i][j];
68 return std::sqrt(x: sum);
69 }
70
71 void setup() {
72
73 N = 3;
74 M1 = M2 = I = Matrix(N,N);
75 M3 = Matrix(3,4);
76 M4 = Matrix(4,3);
77 M5 = Matrix(4, 4, 0.0);
78 M6 = Matrix(4, 4, 0.0);
79
80 M1[0][0] = 1.0; M1[0][1] = 0.9; M1[0][2] = 0.7;
81 M1[1][0] = 0.9; M1[1][1] = 1.0; M1[1][2] = 0.4;
82 M1[2][0] = 0.7; M1[2][1] = 0.4; M1[2][2] = 1.0;
83
84 M2[0][0] = 1.0; M2[0][1] = 0.9; M2[0][2] = 0.7;
85 M2[1][0] = 0.9; M2[1][1] = 1.0; M2[1][2] = 0.3;
86 M2[2][0] = 0.7; M2[2][1] = 0.3; M2[2][2] = 1.0;
87
88 I[0][0] = 1.0; I[0][1] = 0.0; I[0][2] = 0.0;
89 I[1][0] = 0.0; I[1][1] = 1.0; I[1][2] = 0.0;
90 I[2][0] = 0.0; I[2][1] = 0.0; I[2][2] = 1.0;
91
92 M3[0][0] = 1; M3[0][1] = 2; M3[0][2] = 3; M3[0][3] = 4;
93 M3[1][0] = 2; M3[1][1] = 0; M3[1][2] = 2; M3[1][3] = 1;
94 M3[2][0] = 0; M3[2][1] = 1; M3[2][2] = 0; M3[2][3] = 0;
95
96 M4[0][0] = 1; M4[0][1] = 2; M4[0][2] = 400;
97 M4[1][0] = 2; M4[1][1] = 0; M4[1][2] = 1;
98 M4[2][0] = 30; M4[2][1] = 2; M4[2][2] = 0;
99 M4[3][0] = 2; M4[3][1] = 0; M4[3][2] = 1.05;
100
101 // from Higham - nearest correlation matrix
102 M5[0][0] = 2; M5[0][1] = -1; M5[0][2] = 0.0; M5[0][3] = 0.0;
103 M5[1][0] = M5[0][1]; M5[1][1] = 2; M5[1][2] = -1; M5[1][3] = 0.0;
104 M5[2][0] = M5[0][2]; M5[2][1] = M5[1][2]; M5[2][2] = 2; M5[2][3] = -1;
105 M5[3][0] = M5[0][3]; M5[3][1] = M5[1][3]; M5[3][2] = M5[2][3]; M5[3][3] = 2;
106
107 // from Higham - nearest correlation matrix to M5
108 M6[0][0] = 1; M6[0][1] = -0.8084124981; M6[0][2] = 0.1915875019; M6[0][3] = 0.106775049;
109 M6[1][0] = M6[0][1]; M6[1][1] = 1; M6[1][2] = -0.6562326948; M6[1][3] = M6[0][2];
110 M6[2][0] = M6[0][2]; M6[2][1] = M6[1][2]; M6[2][2] = 1; M6[2][3] = M6[0][1];
111 M6[3][0] = M6[0][3]; M6[3][1] = M6[1][3]; M6[3][2] = M6[2][3]; M6[3][3] = 1;
112
113 M7 = M1;
114 M7[0][1] = 0.3; M7[0][2] = 0.2; M7[2][1] = 1.2;
115 }
116
117}
118
119
120void MatricesTest::testEigenvectors() {
121
122 BOOST_TEST_MESSAGE("Testing eigenvalues and eigenvectors calculation...");
123
124 using namespace matrices_test;
125
126 setup();
127
128 Matrix testMatrices[] = { M1, M2 };
129
130 for (auto& M : testMatrices) {
131
132 SymmetricSchurDecomposition dec(M);
133 Array eigenValues = dec.eigenvalues();
134 Matrix eigenVectors = dec.eigenvectors();
135 Real minHolder = QL_MAX_REAL;
136
137 for (Size i=0; i<N; i++) {
138 Array v(N);
139 for (Size j=0; j<N; j++)
140 v[j] = eigenVectors[j][i];
141 // check definition
142 Array a = M*v;
143 Array b = eigenValues[i]*v;
144 if (norm(v: a-b) > 1.0e-15)
145 BOOST_FAIL("Eigenvector definition not satisfied");
146 // check decreasing ordering
147 if (eigenValues[i] >= minHolder) {
148 BOOST_FAIL("Eigenvalues not ordered: " << eigenValues);
149 } else
150 minHolder = eigenValues[i];
151 }
152
153 // check normalization
154 Matrix m = eigenVectors * transpose(m: eigenVectors);
155 if (norm(m: m-I) > 1.0e-15)
156 BOOST_FAIL("Eigenvector not normalized");
157 }
158}
159
160void MatricesTest::testSqrt() {
161
162 BOOST_TEST_MESSAGE("Testing matricial square root...");
163
164 using namespace matrices_test;
165
166 setup();
167
168 Matrix m = pseudoSqrt(M1, SalvagingAlgorithm::None);
169 Matrix temp = m*transpose(m);
170 Real error = norm(m: temp - M1);
171 Real tolerance = 1.0e-12;
172 if (error>tolerance) {
173 BOOST_FAIL("Matrix square root calculation failed\n"
174 << "original matrix:\n" << M1
175 << "pseudoSqrt:\n" << m
176 << "pseudoSqrt*pseudoSqrt:\n" << temp
177 << "\nerror: " << error
178 << "\ntolerance: " << tolerance);
179 }
180}
181
182void MatricesTest::testHighamSqrt() {
183 BOOST_TEST_MESSAGE("Testing Higham matricial square root...");
184
185 using namespace matrices_test;
186
187 setup();
188
189 Matrix tempSqrt = pseudoSqrt(M5, SalvagingAlgorithm::Higham);
190 Matrix ansSqrt = pseudoSqrt(M6, SalvagingAlgorithm::None);
191 Real error = norm(m: ansSqrt - tempSqrt);
192 Real tolerance = 1.0e-4;
193 if (error>tolerance) {
194 BOOST_FAIL("Higham matrix correction failed\n"
195 << "original matrix:\n" << M5
196 << "pseudoSqrt:\n" << tempSqrt
197 << "should be:\n" << ansSqrt
198 << "\nerror: " << error
199 << "\ntolerance: " << tolerance);
200 }
201}
202
203void MatricesTest::testSVD() {
204
205 BOOST_TEST_MESSAGE("Testing singular value decomposition...");
206
207 using namespace matrices_test;
208
209 setup();
210
211 Real tol = 1.0e-12;
212 Matrix testMatrices[] = { M1, M2, M3, M4 };
213
214 for (auto& A : testMatrices) {
215 // m >= n required (rows >= columns)
216 SVD svd(A);
217 // U is m x n
218 const Matrix& U = svd.U();
219 // s is n long
220 Array s = svd.singularValues();
221 // S is n x n
222 Matrix S = svd.S();
223 // V is n x n
224 const Matrix& V = svd.V();
225
226 for (Size i=0; i < S.rows(); i++) {
227 if (S[i][i] != s[i])
228 BOOST_FAIL("S not consistent with s");
229 }
230
231 // tests
232 Matrix U_Utranspose = transpose(m: U)*U;
233 if (norm(m: U_Utranspose-I) > tol)
234 BOOST_FAIL("U not orthogonal (norm of U^T*U-I = "
235 << norm(U_Utranspose-I) << ")");
236
237 Matrix V_Vtranspose = transpose(m: V)*V;
238 if (norm(m: V_Vtranspose-I) > tol)
239 BOOST_FAIL("V not orthogonal (norm of V^T*V-I = "
240 << norm(V_Vtranspose-I) << ")");
241
242 Matrix A_reconstructed = U * S * transpose(m: V);
243 if (norm(m: A_reconstructed-A) > tol)
244 BOOST_FAIL("Product does not recover A: (norm of U*S*V^T-A = "
245 << norm(A_reconstructed-A) << ")");
246 }
247}
248
249void MatricesTest::testQRDecomposition() {
250
251 BOOST_TEST_MESSAGE("Testing QR decomposition...");
252
253 using namespace matrices_test;
254
255 setup();
256
257 Real tol = 1.0e-12;
258 Matrix testMatrices[] = { M1, M2, I,
259 M3, transpose(m: M3), M4, transpose(m: M4), M5 };
260
261 for (const auto& A : testMatrices) {
262 Matrix Q, R;
263 bool pivot = true;
264 const std::vector<Size> ipvt = qrDecomposition(A, q&: Q, r&: R, pivot);
265
266 Matrix P(A.columns(), A.columns(), 0.0);
267
268 // reverse column pivoting
269 for (Size i=0; i < P.columns(); ++i) {
270 P[ipvt[i]][i] = 1.0;
271 }
272
273 if (norm(m: Q*R - A*P) > tol)
274 BOOST_FAIL("Q*R does not match matrix A*P (norm = "
275 << norm(Q*R-A*P) << ")");
276
277 pivot = false;
278 qrDecomposition(A, q&: Q, r&: R, pivot);
279
280 if (norm(m: Q*R - A) > tol)
281 BOOST_FAIL("Q*R does not match matrix A (norm = "
282 << norm(Q*R-A) << ")");
283 }
284}
285
286void MatricesTest::testQRSolve() {
287
288 BOOST_TEST_MESSAGE("Testing QR solve...");
289
290 using namespace matrices_test;
291
292 setup();
293
294 Real tol = 1.0e-12;
295 MersenneTwisterUniformRng rng(1234);
296 Matrix bigM(50, 100, 0.0);
297 for (Size i=0; i < std::min(a: bigM.rows(), b: bigM.columns()); ++i) {
298 bigM[i][i] = i+1.0;
299 }
300
301 Matrix randM(50, 200);
302 for (Size i=0; i < randM.rows(); ++i)
303 for (Size j=0; j < randM.columns(); ++j)
304 randM[i][j] = rng.next().value;
305
306 Matrix testMatrices[] = {M1, M2, M3, transpose(m: M3),
307 M4, transpose(m: M4), M5, I, M7,
308 bigM, transpose(m: bigM),
309 randM, transpose(m: randM) };
310
311 for (const auto& A : testMatrices) {
312 Array b(A.rows());
313
314 for (Size k=0; k < 10; ++k) {
315 for (Real& iter : b) {
316 iter = rng.next().value;
317 }
318 const Array x = qrSolve(a: A, b, pivot: true);
319
320 if (A.columns() >= A.rows()) {
321 if (norm(v: A*x - b) > tol)
322 BOOST_FAIL("A*x does not match vector b (norm = "
323 << norm(A*x - b) << ")");
324 }
325 else {
326 // use the SVD to calculate the reference values
327 const Size n = A.columns();
328 Array xr(n, 0.0);
329
330 SVD svd(A);
331 const Matrix& V = svd.V();
332 const Matrix& U = svd.U();
333 const Array& w = svd.singularValues();
334 const Real threshold = n*QL_EPSILON;
335
336 for (Size i=0; i<n; ++i) {
337 if (w[i] > threshold) {
338 const Real u = std::inner_product(first1: U.column_begin(i),
339 last1: U.column_end(i),
340 first2: b.begin(), init: Real(0.0))/w[i];
341
342 for (Size j=0; j<n; ++j) {
343 xr[j] +=u*V[j][i];
344 }
345 }
346 }
347
348 if (norm(v: xr-x) > tol) {
349 BOOST_FAIL("least square solution does not match (norm = "
350 << norm(x - xr) << ")");
351
352 }
353 }
354 }
355 }
356}
357
358void MatricesTest::testInverse() {
359
360 BOOST_TEST_MESSAGE("Testing LU inverse calculation...");
361
362 using namespace matrices_test;
363
364 setup();
365
366 Real tol = 1.0e-12;
367 Matrix testMatrices[] = { M1, M2, I, M5 };
368
369 for (const auto& A : testMatrices) {
370 const Matrix invA = inverse(m: A);
371
372 const Matrix I1 = invA*A;
373 const Matrix I2 = A*invA;
374
375 Matrix identity(A.rows(), A.rows(), 0.0);
376 for (Size i=0; i < A.rows(); ++i) identity[i][i] = 1.0;
377
378 if (norm(m: I1 - identity) > tol)
379 BOOST_FAIL("inverse(A)*A does not recover unit matrix (norm = "
380 << norm(I1-identity) << ")");
381
382 if (norm(m: I2 - identity) > tol)
383 BOOST_FAIL("A*inverse(A) does not recover unit matrix (norm = "
384 << norm(I1-identity) << ")");
385 }
386}
387
388void MatricesTest::testDeterminant() {
389
390 BOOST_TEST_MESSAGE("Testing LU determinant calculation...");
391
392 using namespace matrices_test;
393
394 setup();
395 Real tol = 1e-10;
396
397 Matrix testMatrices[] = {M1, M2, M5, M6, I};
398 // expected results calculated with octave
399 Real expected[] = { 0.044, -0.012, 5.0, 5.7621e-11, 1.0};
400
401 for (Size j=0; j<LENGTH(testMatrices); ++j) {
402 const Real calculated = determinant(m: testMatrices[j]);
403 if (std::fabs(x: expected[j] - calculated) > tol)
404 BOOST_FAIL("determinant calculation failed "
405 << "\n matrix :\n" << testMatrices[j]
406 << "\n calculated : " << calculated
407 << "\n expected : " << expected[j]);
408 }
409
410 MersenneTwisterUniformRng rng(1234);
411 for (Size j=0; j<100; ++j) {
412 Matrix m(3, 3, 0.0);
413 for (Real& iter : m)
414 iter = rng.next().value;
415
416 if ((j % 3) == 0U) {
417 // every third matrix is a singular matrix
418 Size row = Size(3*rng.next().value);
419 std::fill(first: m.row_begin(i: row), last: m.row_end(i: row), value: 0.0);
420 }
421
422 Real a=m[0][0];
423 Real b=m[0][1];
424 Real c=m[0][2];
425 Real d=m[1][0];
426 Real e=m[1][1];
427 Real f=m[1][2];
428 Real g=m[2][0];
429 Real h=m[2][1];
430 Real i=m[2][2];
431
432 const Real expected = a*e*i+b*f*g+c*d*h-(g*e*c+h*f*a+i*d*b);
433 const Real calculated = determinant(m);
434
435 if (std::fabs(x: expected-calculated) > tol)
436 BOOST_FAIL("determinant calculation failed "
437 << "\n matrix :\n" << m
438 << "\n calculated : " << calculated
439 << "\n expected : " << expected);
440 }
441}
442
443void MatricesTest::testOrthogonalProjection() {
444 BOOST_TEST_MESSAGE("Testing orthogonal projections...");
445
446 Size dimension = 1000;
447 Size numberVectors = 50;
448 Real multiplier = 100;
449 Real tolerance = 1e-6;
450 unsigned long seed = 1;
451
452 Real errorAcceptable = 1E-11;
453
454 Matrix test(numberVectors,dimension);
455
456 MersenneTwisterUniformRng rng(seed);
457
458 for (Size i=0; i < numberVectors; ++i)
459 for (Size j=0; j < dimension; ++j)
460 test[i][j] = rng.next().value;
461
462 OrthogonalProjections projector(test,
463 multiplier,
464 tolerance );
465
466 Size numberFailures =0;
467 Size failuresTwo=0;
468
469 for (Size i=0; i < numberVectors; ++i)
470 {
471 // check output vector i is orthogonal to all other vectors
472
473 if (projector.validVectors()[i])
474 {
475 for (Size j=0; j < numberVectors; ++j)
476 if (projector.validVectors()[j] && i != j)
477 {
478 Real dotProduct=0.0;
479 for (Size k=0; k < dimension; ++k)
480 dotProduct += test[j][k]*projector.GetVector(index: i)[k];
481
482 if (fabs(x: dotProduct) > errorAcceptable)
483 ++numberFailures;
484
485 }
486
487 Real innerProductWithOriginal =0.0;
488 Real normSq =0.0;
489
490 for (Size j=0; j < dimension; ++j)
491 {
492 innerProductWithOriginal += projector.GetVector(index: i)[j]*test[i][j];
493 normSq += test[i][j]*test[i][j];
494 }
495
496 if (fabs(x: innerProductWithOriginal-normSq) > errorAcceptable)
497 ++failuresTwo;
498
499 }
500
501 }
502
503 if (numberFailures > 0 || failuresTwo >0)
504 BOOST_FAIL("OrthogonalProjections test failed with " << numberFailures << " failures of orthogonality and "
505 << failuresTwo << " failures of projection size.");
506
507}
508
509void MatricesTest::testCholeskyDecomposition() {
510
511 BOOST_TEST_MESSAGE("Testing Cholesky Decomposition...");
512
513 // This test case fails prior to release 1.8
514
515 // The eigenvalues of this matrix are
516 // 0.0438523; 0.0187376; 0.000245617; 0.000127656; 8.35899e-05; 6.14215e-05;
517 // 1.94241e-05; 1.14417e-06; 9.79481e-18; 1.31141e-18; 5.81155e-19
518
519 Real tmp[11][11] = {
520 {6.4e-05, 5.28e-05, 2.28e-05, 0.00032, 0.00036, 6.4e-05,
521 6.3968010664e-06, 7.2e-05, 7.19460269899e-06, 1.2e-05,
522 1.19970004999e-06},
523 {5.28e-05, 0.000121, 1.045e-05, 0.00044, 0.000165, 2.2e-05,
524 2.19890036657e-06, 1.65e-05, 1.64876311852e-06, 1.1e-05,
525 1.09972504583e-06},
526 {2.28e-05, 1.045e-05, 9.025e-05, 0, 0.0001425, 9.5e-06,
527 9.49525158294e-07, 2.85e-05, 2.84786356835e-06, 4.75e-06,
528 4.74881269789e-07},
529 {0.00032, 0.00044, 0, 0.04, 0.009, 0.0008, 7.996001333e-05, 0.0006,
530 5.99550224916e-05, 0.0001, 9.99750041661e-06},
531 {0.00036, 0.000165, 0.0001425, 0.009, 0.0225, 0.0003, 2.99850049987e-05,
532 0.001125, 0.000112415667172, 0.000225, 2.24943759374e-05},
533 {6.4e-05, 2.2e-05, 9.5e-06, 0.0008, 0.0003, 0.0001, 9.99500166625e-06,
534 7.5e-05, 7.49437781145e-06, 2e-05, 1.99950008332e-06},
535 {6.3968010664e-06, 2.19890036657e-06, 9.49525158294e-07,
536 7.996001333e-05, 2.99850049987e-05, 9.99500166625e-06,
537 9.99000583083e-07, 7.49625124969e-06, 7.49063187129e-07,
538 1.99900033325e-06, 1.99850066645e-07},
539 {7.2e-05, 1.65e-05, 2.85e-05, 0.0006, 0.001125, 7.5e-05,
540 7.49625124969e-06, 0.000225, 2.24831334343e-05, 1.5e-05,
541 1.49962506249e-06},
542 {7.19460269899e-06, 1.64876311852e-06, 2.84786356835e-06,
543 5.99550224916e-05, 0.000112415667172, 7.49437781145e-06,
544 7.49063187129e-07, 2.24831334343e-05, 2.24662795123e-06,
545 1.49887556229e-06, 1.49850090584e-07},
546 {1.2e-05, 1.1e-05, 4.75e-06, 0.0001, 0.000225, 2e-05, 1.99900033325e-06,
547 1.5e-05, 1.49887556229e-06, 2.5e-05, 2.49937510415e-06},
548 {1.19970004999e-06, 1.09972504583e-06, 4.74881269789e-07,
549 9.99750041661e-06, 2.24943759374e-05, 1.99950008332e-06,
550 1.99850066645e-07, 1.49962506249e-06, 1.49850090584e-07,
551 2.49937510415e-06, 2.49875036451e-07}};
552
553 Matrix m(11,11);
554 for(Size i=0;i<11;++i) {
555 for(Size j=0;j<11;++j) {
556 m[i][j] = tmp[i][j];
557 }
558 }
559
560 Matrix c = CholeskyDecomposition(m,flexible: true);
561 Matrix m2 = c * transpose(m: c);
562
563 Real tol = 1.0E-12;
564 for(Size i=0;i<11;++i) {
565 for(Size j=0;j<11;++j) {
566 if(std::isnan(x: m2[i][j])) {
567 BOOST_FAIL("Faield to verify Cholesky decomposition at (i,j)=("
568 << i << "," << j << "), replicated value is nan");
569 }
570 // this does not detect nan values
571 if(std::abs(x: m[i][j]-m2[i][j]) > tol) {
572 BOOST_FAIL("Failed to verify Cholesky decomposition at (i,j)=("
573 << i << "," << j << "), original value is "
574 << m[i][j] << ", replicated value is " << m2[i][j]);
575 }
576 }
577 }
578}
579
580void MatricesTest::testMoorePenroseInverse() {
581
582 BOOST_TEST_MESSAGE("Testing Moore-Penrose inverse...");
583
584 // this is taken from
585 // http://de.mathworks.com/help/matlab/ref/pinv.html
586 Real tmp[8][6] = {{64, 2, 3, 61, 60, 6}, {9, 55, 54, 12, 13, 51},
587 {17, 47, 46, 20, 21, 43}, {40, 26, 27, 37, 36, 30},
588 {32, 34, 35, 29, 28, 38}, {41, 23, 22, 44, 45, 19},
589 {49, 15, 14, 52, 53, 11}, {8, 58, 59, 5, 4, 62}};
590 Matrix A(8, 6);
591 for (Size i = 0; i < 8; ++i) {
592 for (Size j = 0; j < 6; ++j) {
593 A(i, j) = tmp[i][j];
594 }
595 }
596
597 Matrix P = moorePenroseInverse(A);
598 Array b(8, 260.0);
599 Array x = P*b;
600
601 Real cached[6] = {1.153846153846152, 1.461538461538463, 1.384615384615384,
602 1.384615384615385, 1.461538461538462, 1.153846153846152};
603 constexpr double tol = 500.0 * QL_EPSILON;
604
605 for (Size i = 0; i < 6; ++i) {
606 if (std::abs(x: x[i] - cached[i]) > tol) {
607 BOOST_FAIL("Failed to verify minimal norm solution obtained from "
608 "Moore-Penrose-Inverse against cached results, component "
609 << i << " is " << x[i] << ", expected " << cached[i]
610 << ", difference " << x[i] - cached[i] << ", tolerance "
611 << tol);
612 }
613 }
614
615 Array y = A*x;
616 constexpr double tol2 = 2000.0 * QL_EPSILON;
617 for (Size i = 0; i < 6; ++i) {
618 if (std::abs(x: y[i] - 260.0) > tol2) {
619 BOOST_FAIL(
620 "Failed to verify minimal norm solution obtained from "
621 "Moore-Penrose-Inverse when back-substituting, rhs component "
622 << i << " is " << y[i] << ", expected 260.0, difference "
623 << y[i] - 260.0 << ", tolerance " << tol2);
624 }
625 }
626
627}
628
629
630namespace matrices_test {
631 class MatrixMult {
632 public:
633 explicit MatrixMult(Matrix m) : m_(std::move(m)) {}
634 Array operator()(const Array& x) const {
635 return m_ * x;
636 }
637
638 private:
639 const Matrix m_;
640 };
641
642 Real norm2(const Array& x) {
643 return std::sqrt(x: DotProduct(v1: x,v2: x));
644 }
645}
646
647void MatricesTest::testIterativeSolvers() {
648 BOOST_TEST_MESSAGE("Testing iterative solvers...");
649
650 using namespace matrices_test;
651
652 setup();
653
654 Array b(3);
655 b[0] = 1.0; b[1] = 0.5; b[2] = 3.0;
656
657 constexpr double relTol = 1e4 * QL_EPSILON;
658
659 const Array x = BiCGstab(MatrixMult(M1), 3, relTol).solve(b).x;
660 if (norm2(x: M1*x-b)/norm2(x: b) > relTol) {
661 BOOST_FAIL("Failed to calculate inverse using BiCGstab"
662 << "\n rel error : " << norm2(M1*x-b)/norm2(b)
663 << "\n rel tolerance : " << relTol);
664 }
665
666 const GMRESResult u = GMRES(MatrixMult(M1), 3, relTol).solve(b, x0: b);
667 if (norm2(x: M1*u.x-b)/norm2(x: b) > relTol) {
668 BOOST_FAIL("Failed to calculate inverse using gmres"
669 << "\n rel error : " << norm2(M1*u.x-b)/norm2(b)
670 << "\n rel tolerance : " << relTol);
671 }
672 const Array errors = Array(u.errors.begin(), u.errors.end());
673 for (Real error : errors) {
674 const Array x = GMRES(MatrixMult(M1), 10, 1.01 * error).solve(b, x0: b).x;
675
676 const Real calculated = norm2(x: M1*x-b)/norm2(x: b);
677 const Real expected = error;
678
679 if (std::fabs(x: calculated - expected) > relTol) {
680 BOOST_FAIL("Failed to calculate solution error"
681 << "\n calculated error: " << calculated
682 << "\n expected error : " << expected);
683 }
684 }
685
686 const Array v = GMRES(MatrixMult(M1), 1, relTol,
687 MatrixMult(inverse(m: M1))).solve(b, x0: b).x;
688
689 if (norm2(x: M1*v-b)/norm2(x: b) > relTol) {
690 BOOST_FAIL("Failed to calculate inverse using gmres "
691 "with exact preconditioning"
692 << "\n rel error : " << norm2(M1*v-b)/norm2(b)
693 << "\n rel tolerance : " << relTol);
694 }
695
696 const Array w = GMRES(MatrixMult(M1), 3, relTol,
697 MatrixMult(M1)).solve(b, x0: b).x;
698 if (norm2(x: M1*w-b)/norm2(x: b) > relTol) {
699 BOOST_FAIL("Failed to calculate inverse using gmres "
700 "with nonsense preconditioning"
701 << "\n rel error : " << norm2(M1*w-b)/norm2(b)
702 << "\n rel tolerance : " << relTol);
703 }
704}
705
706void MatricesTest::testInitializers() {
707 BOOST_TEST_MESSAGE("Testing matrix initializers...");
708
709 Matrix m1 = {};
710 BOOST_REQUIRE(m1.rows() == 0);
711 BOOST_REQUIRE(m1.columns() == 0);
712
713 Matrix m2 = {
714 {1.0, 2.0, 3.0},
715 {4.0, 5.0, 6.0}
716 };
717 BOOST_REQUIRE(m2.rows() == 2);
718 BOOST_REQUIRE(m2.columns() == 3);
719 BOOST_CHECK_EQUAL(m2(0, 0), 1.0);
720 BOOST_CHECK_EQUAL(m2(0, 1), 2.0);
721 BOOST_CHECK_EQUAL(m2(0, 2), 3.0);
722 BOOST_CHECK_EQUAL(m2(1, 0), 4.0);
723 BOOST_CHECK_EQUAL(m2(1, 1), 5.0);
724 BOOST_CHECK_EQUAL(m2(1, 2), 6.0);
725}
726
727
728namespace {
729
730 typedef std::pair< std::pair< std::vector<Size>, std::vector<Size> >,
731 std::vector<Real> > coordinate_tuple;
732
733 coordinate_tuple sparseMatrixToCoordinateTuple(const SparseMatrix& m) {
734 std::vector<Size> row_idx, col_idx;
735 std::vector<Real> data;
736 for (auto iter1 = m.begin1(); iter1 != m.end1(); ++iter1)
737 for (auto iter2 = iter1.begin(); iter2 != iter1.end(); ++iter2) {
738 row_idx.push_back(x: iter1.index1());
739 col_idx.push_back(x: iter2.index2());
740 data.push_back(x: *iter2);
741 }
742
743 return std::make_pair(x: std::make_pair(x&: row_idx, y&: col_idx), y&: data);
744 }
745
746}
747
748void MatricesTest::testSparseMatrixMemory() {
749
750 BOOST_TEST_MESSAGE("Testing sparse matrix memory layout...");
751
752 SparseMatrix m(8, 4);
753 BOOST_CHECK_EQUAL(m.filled1(), 1);
754 BOOST_CHECK_EQUAL(m.size1(), 8);
755 BOOST_CHECK_EQUAL(m.size2(), 4);
756 BOOST_CHECK_EQUAL(std::distance(m.begin1(), m.end1()), m.size1());
757
758 auto coords = sparseMatrixToCoordinateTuple(m);
759 BOOST_CHECK_EQUAL(coords.first.first.size(), 0);
760
761 m(3, 1) = 42;
762 coords = sparseMatrixToCoordinateTuple(m);
763 BOOST_CHECK_EQUAL(std::distance(m.begin1(), m.end1()), m.size1());
764 BOOST_CHECK_EQUAL(coords.first.first.size(), 1);
765 BOOST_CHECK_EQUAL(coords.first.first[0], 3);
766 BOOST_CHECK_EQUAL(coords.first.second[0], 1);
767 BOOST_CHECK_EQUAL(coords.second[0], 42);
768
769 m(1, 2) = 6;
770 coords = sparseMatrixToCoordinateTuple(m);
771 BOOST_CHECK_EQUAL(coords.first.first.size(), 2);
772 BOOST_CHECK_EQUAL(coords.first.first[0], 1);
773 BOOST_CHECK_EQUAL(coords.first.second[0], 2);
774 BOOST_CHECK_EQUAL(coords.second[0], 6);
775
776 Array x{1, 2, 3, 4};
777 Array y = prod(A: m, x);
778 BOOST_CHECK_EQUAL(y, Array({0, 18, 0, 84}));
779
780 m(3, 2) = 43;
781 coords = sparseMatrixToCoordinateTuple(m);
782 BOOST_CHECK_EQUAL(coords.first.first.size(), 3);
783 BOOST_CHECK_EQUAL(coords.first.first[2], 3);
784 BOOST_CHECK_EQUAL(coords.first.second[2], 2);
785 BOOST_CHECK_EQUAL(coords.second[2], 43);
786
787 m(7, 3) = 44;
788 coords = sparseMatrixToCoordinateTuple(m);
789 BOOST_CHECK_EQUAL(coords.first.first.size(), 4);
790 BOOST_CHECK_EQUAL(coords.first.first[3], 7);
791 BOOST_CHECK_EQUAL(coords.first.second[3], 3);
792 BOOST_CHECK_EQUAL(coords.second[3], 44);
793
794 Size entries(0);
795 for (auto iter1 = m.begin1(); iter1 != m.end1(); ++iter1)
796 entries+=std::distance(first: iter1.begin(), last: iter1.end());
797
798 BOOST_CHECK_EQUAL(entries, 4);
799
800}
801
802#define QL_CHECK_CLOSE_MATRIX(actual, expected) \
803 BOOST_REQUIRE(actual.rows() == expected.rows() && \
804 actual.columns() == expected.columns()); \
805 for (auto i = 0u; i < actual.rows(); i++) { \
806 for (auto j = 0u; j < actual.columns(); j++) { \
807 QL_CHECK_CLOSE(actual(i, j), expected(i, j), 100 * QL_EPSILON); \
808 } \
809 } \
810
811void MatricesTest::testOperators() {
812
813 BOOST_TEST_MESSAGE("Testing matrix operators...");
814
815 auto get_matrix = []() {
816 return Matrix(2, 3, 4.0);
817 };
818
819 const auto m = get_matrix();
820
821 const auto negative = Matrix(2, 3, -4.0);
822 const auto lvalue_negative = -m;
823 const auto rvalue_negative = -get_matrix();
824
825 QL_CHECK_CLOSE_MATRIX(lvalue_negative, negative);
826 QL_CHECK_CLOSE_MATRIX(rvalue_negative, negative);
827
828 const auto matrix_sum = Matrix(2, 3, 8.0);
829 const auto lvalue_lvalue_sum = m + m;
830 const auto lvalue_rvalue_sum = m + get_matrix();
831 const auto rvalue_lvalue_sum = get_matrix() + m;
832 const auto rvalue_rvalue_sum = get_matrix() + get_matrix();
833
834 QL_CHECK_CLOSE_MATRIX(lvalue_lvalue_sum, matrix_sum);
835 QL_CHECK_CLOSE_MATRIX(lvalue_rvalue_sum, matrix_sum);
836 QL_CHECK_CLOSE_MATRIX(rvalue_lvalue_sum, matrix_sum);
837 QL_CHECK_CLOSE_MATRIX(rvalue_rvalue_sum, matrix_sum);
838
839 const auto matrix_difference = Matrix(2, 3, 0.0);
840 const auto lvalue_lvalue_difference = m - m; // NOLINT(misc-redundant-expression)
841 const auto lvalue_rvalue_difference = m - get_matrix();
842 const auto rvalue_lvalue_difference = get_matrix() - m;
843 const auto rvalue_rvalue_difference = get_matrix() - get_matrix();
844
845 QL_CHECK_CLOSE_MATRIX(lvalue_lvalue_difference, matrix_difference);
846 QL_CHECK_CLOSE_MATRIX(lvalue_rvalue_difference, matrix_difference);
847 QL_CHECK_CLOSE_MATRIX(rvalue_lvalue_difference, matrix_difference);
848 QL_CHECK_CLOSE_MATRIX(rvalue_rvalue_difference, matrix_difference);
849
850 const auto scalar_product = Matrix(2, 3, 6.0);
851 const auto lvalue_real_product = m * 1.5;
852 const auto rvalue_real_product = get_matrix() * 1.5;
853 const auto real_lvalue_product = 1.5 * m;
854 const auto real_rvalue_product = 1.5 * get_matrix();
855
856 QL_CHECK_CLOSE_MATRIX(lvalue_real_product, scalar_product);
857 QL_CHECK_CLOSE_MATRIX(rvalue_real_product, scalar_product);
858 QL_CHECK_CLOSE_MATRIX(real_lvalue_product, scalar_product);
859 QL_CHECK_CLOSE_MATRIX(real_rvalue_product, scalar_product);
860
861 const auto scalar_quotient = Matrix(2, 3, 2.0);
862 const auto lvalue_real_quotient = m / 2.0;
863 const auto rvalue_real_quotient = get_matrix() / 2.0;
864
865 QL_CHECK_CLOSE_MATRIX(lvalue_real_quotient, scalar_quotient);
866 QL_CHECK_CLOSE_MATRIX(rvalue_real_quotient, scalar_quotient);
867}
868
869test_suite* MatricesTest::suite() {
870 auto* suite = BOOST_TEST_SUITE("Matrix tests");
871
872 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testOrthogonalProjection));
873 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testEigenvectors));
874 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testSqrt));
875 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testSVD));
876 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testHighamSqrt));
877 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testQRDecomposition));
878 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testQRSolve));
879 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testInverse));
880 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testDeterminant));
881 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testSparseMatrixMemory));
882 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testCholeskyDecomposition));
883 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testMoorePenroseInverse));
884 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testIterativeSolvers));
885 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testInitializers));
886 suite->add(QUANTLIB_TEST_CASE(&MatricesTest::testOperators));
887 return suite;
888}
889
890

source code of quantlib/test-suite/matrices.cpp