test/inverse_svd.cc

Code
Comments
Other
Rev Date Author Line
4282 29 Jan 23 peter 1 // $Id$
4282 29 Jan 23 peter 2
4282 29 Jan 23 peter 3 /*
4282 29 Jan 23 peter 4   Copyright (C) 2023 Peter Johansson
4282 29 Jan 23 peter 5
4282 29 Jan 23 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4282 29 Jan 23 peter 7
4282 29 Jan 23 peter 8   The yat library is free software; you can redistribute it and/or
4282 29 Jan 23 peter 9   modify it under the terms of the GNU General Public License as
4282 29 Jan 23 peter 10   published by the Free Software Foundation; either version 3 of the
4282 29 Jan 23 peter 11   License, or (at your option) any later version.
4282 29 Jan 23 peter 12
4282 29 Jan 23 peter 13   The yat library is distributed in the hope that it will be useful,
4282 29 Jan 23 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4282 29 Jan 23 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4282 29 Jan 23 peter 16   General Public License for more details.
4282 29 Jan 23 peter 17
4282 29 Jan 23 peter 18   You should have received a copy of the GNU General Public License
4282 29 Jan 23 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4282 29 Jan 23 peter 20 */
4282 29 Jan 23 peter 21
4282 29 Jan 23 peter 22 #include <config.h>
4282 29 Jan 23 peter 23
4282 29 Jan 23 peter 24 #include "Suite.h"
4282 29 Jan 23 peter 25
4283 30 Jan 23 peter 26 #include "yat/utility/DiagonalMatrix.h"
4282 29 Jan 23 peter 27 #include "yat/utility/Matrix.h"
4283 30 Jan 23 peter 28 #include "yat/utility/SVD.h"
4282 29 Jan 23 peter 29
4378 11 Oct 23 peter 30 #include <cmath>
4378 11 Oct 23 peter 31
4283 30 Jan 23 peter 32 using namespace theplu::yat;
4283 30 Jan 23 peter 33 using utility::Matrix;
4283 30 Jan 23 peter 34
4283 30 Jan 23 peter 35 bool equal(test::Suite& suite, const Matrix& A, const Matrix& B, double margin)
4283 30 Jan 23 peter 36 {
4283 30 Jan 23 peter 37   if (A.rows() != B.rows() || A.columns() != B.columns()) {
4283 30 Jan 23 peter 38     suite.add(false);
4283 30 Jan 23 peter 39     suite.err() << "error: incorrect dimensions: "
4283 30 Jan 23 peter 40                 << A.rows() << " x " << A.columns() << " vs "
4283 30 Jan 23 peter 41                 << B.rows() << " x " << B.columns() << "\n";
4283 30 Jan 23 peter 42     return false;
4283 30 Jan 23 peter 43   }
4283 30 Jan 23 peter 44   bool ok = true;
4283 30 Jan 23 peter 45   for (size_t i=0; i<A.rows(); ++i)
4283 30 Jan 23 peter 46     for (size_t j=0; j<A.columns(); ++j)
4283 30 Jan 23 peter 47       if (!suite.equal_fix(A(i,j), B(i,j), margin)) {
4283 30 Jan 23 peter 48         suite.err() << "error: element: " << i << ", " << j << "\n";
4283 30 Jan 23 peter 49         ok = false;
4283 30 Jan 23 peter 50       }
4283 30 Jan 23 peter 51   suite.add(ok);
4283 30 Jan 23 peter 52   return ok;
4283 30 Jan 23 peter 53 }
4283 30 Jan 23 peter 54
4283 30 Jan 23 peter 55
4282 29 Jan 23 peter 56 int main(int argc, char* argv[])
4282 29 Jan 23 peter 57 {
4282 29 Jan 23 peter 58   test::Suite suite(argc, argv);
4282 29 Jan 23 peter 59
4282 29 Jan 23 peter 60   size_t n = 4;
4283 30 Jan 23 peter 61   Matrix Eye(n, n, 0.0);
4282 29 Jan 23 peter 62   for (size_t i=0; i<n; ++i)
4283 30 Jan 23 peter 63     Eye(i,i) = 1.0;
4283 30 Jan 23 peter 64
4285 30 Jan 23 peter 65   Matrix A = test::generate_Matrix(n, n, 0.0);
4283 30 Jan 23 peter 66   // avoid near-zero eigenvalues
4283 30 Jan 23 peter 67   A += 0.1 * Eye;
4283 30 Jan 23 peter 68   Matrix B;
4282 29 Jan 23 peter 69   inverse_svd(A, B);
4283 30 Jan 23 peter 70
4283 30 Jan 23 peter 71
4282 29 Jan 23 peter 72   if (B.rows() != n || B.columns() != n) {
4282 29 Jan 23 peter 73     suite.add(false);
4282 29 Jan 23 peter 74     suite.err() << "inverse_svd: incorrect dimensions: "
4282 29 Jan 23 peter 75                 << B.rows() << " x " << B.columns() << "\n";
4282 29 Jan 23 peter 76   }
4282 29 Jan 23 peter 77   else {
4283 30 Jan 23 peter 78
4283 30 Jan 23 peter 79     utility::SVD svd(A);
4283 30 Jan 23 peter 80     svd.decompose();
4283 30 Jan 23 peter 81     utility::DiagonalMatrix D(svd.s());
4283 30 Jan 23 peter 82     Matrix UDVt = svd.U() * D * transpose(svd.V());
4283 30 Jan 23 peter 83     suite.out() << "compare A and UDVt\n";
4283 30 Jan 23 peter 84     equal(suite, A, UDVt, 1e-9);
4283 30 Jan 23 peter 85
4283 30 Jan 23 peter 86     const Matrix VtV = transpose(svd.V()) * svd.V();
4283 30 Jan 23 peter 87     suite.out() << "compare VtV and Eye\n";
4283 30 Jan 23 peter 88     equal(suite, VtV, Eye, 1e-9);
4283 30 Jan 23 peter 89
4283 30 Jan 23 peter 90     Matrix UUt = svd.U() * transpose(svd.U());
4283 30 Jan 23 peter 91     suite.out() << "compare UUt and Eye\n";
4283 30 Jan 23 peter 92     equal(suite, UUt, Eye, 1e-9);
4283 30 Jan 23 peter 93
4283 30 Jan 23 peter 94     utility::DiagonalMatrix D_inv(n, n);
4282 29 Jan 23 peter 95     for (size_t i=0; i<n; ++i)
4283 30 Jan 23 peter 96       D_inv(i) = 1.0 / D(i,i);
4282 29 Jan 23 peter 97
4283 30 Jan 23 peter 98     Matrix X1 = (D * Eye) * D_inv;
4283 30 Jan 23 peter 99     suite.out() << "compare X1 and Eye\n";
4283 30 Jan 23 peter 100     equal(suite, X1, Eye, 1e-9);
4283 30 Jan 23 peter 101
4283 30 Jan 23 peter 102     Matrix X2 = (D * VtV) * D_inv;
4283 30 Jan 23 peter 103     suite.out() << "compare X2 and Eye\n";
4283 30 Jan 23 peter 104     equal(suite, X2, Eye, 1e-9);
4283 30 Jan 23 peter 105
4283 30 Jan 23 peter 106     Matrix X3 = D * (VtV * D_inv);
4283 30 Jan 23 peter 107     suite.out() << "compare X3 and Eye\n";
4283 30 Jan 23 peter 108     equal(suite, X3, Eye, 1e-9);
4283 30 Jan 23 peter 109
4283 30 Jan 23 peter 110
4283 30 Jan 23 peter 111     suite.out() << "test that B is inverse of A\n";
4283 30 Jan 23 peter 112     Matrix C = A * B;
4283 30 Jan 23 peter 113
4283 30 Jan 23 peter 114     suite.out() << "compare C and Eye\n";
4283 30 Jan 23 peter 115     equal(suite, C, Eye, 1e-9);
4282 29 Jan 23 peter 116   }
4282 29 Jan 23 peter 117
4282 29 Jan 23 peter 118   return suite.return_value();
4282 29 Jan 23 peter 119 }

Generated on Sat Oct 14 02:21:01 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)