test/svd.cc

Code
Comments
Other
Rev Date Author Line
4 25 Feb 03 daniel 1 // $Id$
4 25 Feb 03 daniel 2
675 10 Oct 06 jari 3 /*
831 27 Mar 07 peter 4   Copyright (C) 2003 Daniel Dalevi
2119 12 Dec 09 peter 5   Copyright (C) 2004 Jari Häkkinen
2119 12 Dec 09 peter 6   Copyright (C) 2005 Jari Häkkinen, Peter Johansson
2119 12 Dec 09 peter 7   Copyright (C) 2006 Jari Häkkinen
2119 12 Dec 09 peter 8   Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 9   Copyright (C) 2012 Peter Johansson
4 25 Feb 03 daniel 10
1437 25 Aug 08 peter 11   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 12
675 10 Oct 06 jari 13   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 14   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 15   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 16   License, or (at your option) any later version.
675 10 Oct 06 jari 17
675 10 Oct 06 jari 18   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 19   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 20   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 21   General Public License for more details.
675 10 Oct 06 jari 22
675 10 Oct 06 jari 23   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 24   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 25 */
675 10 Oct 06 jari 26
2881 18 Nov 12 peter 27 #include <config.h>
2881 18 Nov 12 peter 28
1248 19 Mar 08 peter 29 #include "Suite.h"
1248 19 Mar 08 peter 30
675 10 Oct 06 jari 31 #include "yat/random/random.h"
1121 22 Feb 08 peter 32 #include "yat/utility/Matrix.h"
675 10 Oct 06 jari 33 #include "yat/utility/SVD.h"
1120 21 Feb 08 peter 34 #include "yat/utility/Vector.h"
675 10 Oct 06 jari 35
680 11 Oct 06 jari 36 using namespace theplu::yat;
4 25 Feb 03 daniel 37
1248 19 Mar 08 peter 38 void do_test(size_t m, size_t n, utility::SVD::SVDalgorithm algo,
1248 19 Mar 08 peter 39              test::Suite& suite)
4 25 Feb 03 daniel 40 {
1248 19 Mar 08 peter 41   // initialise a random test-matrix
1248 19 Mar 08 peter 42   theplu::yat::random::ContinuousUniform rnd;
227 01 Feb 05 jari 43
1248 19 Mar 08 peter 44   // set seed to make test deterministic
1248 19 Mar 08 peter 45   theplu::yat::random::RNG::instance()->seed(20);
227 01 Feb 05 jari 46
1121 22 Feb 08 peter 47   utility::Matrix A(m,n);
227 01 Feb 05 jari 48   for (size_t i=0; i<m; ++i)
227 01 Feb 05 jari 49     for(size_t j=0; j<n; ++j)
377 07 Aug 05 jari 50       A(i,j)=1000*rnd();
227 01 Feb 05 jari 51
301 30 Apr 05 peter 52   utility::SVD svd(A);
227 01 Feb 05 jari 53    svd.decompose(algo);
1120 21 Feb 08 peter 54   theplu::yat::utility::Vector s(svd.s());
1121 22 Feb 08 peter 55   utility::Matrix S(s.size(),s.size());
1476 04 Sep 08 peter 56   for (size_t i=1; i<s.size(); ++i)
1476 04 Sep 08 peter 57     if (!suite.add(s(i-1)>=s(i)))
1476 04 Sep 08 peter 58       suite.err() << "single values not in non-increasing order\n";
227 01 Feb 05 jari 59   for (size_t i=0; i<s.size(); ++i)
1151 25 Feb 08 peter 60     S(i,i)=s(i);
1121 22 Feb 08 peter 61   utility::Matrix Vtranspose=svd.V();
420 02 Dec 05 jari 62   Vtranspose.transpose();
420 02 Dec 05 jari 63   // Reconstructing A = U*S*Vtranspose
1121 22 Feb 08 peter 64   utility::Matrix Areconstruct(svd.U());
420 02 Dec 05 jari 65   Areconstruct*=S;
420 02 Dec 05 jari 66   Areconstruct*=Vtranspose;
1251 03 Apr 08 peter 67   if (!suite.equal_range(A.begin(), A.end(), Areconstruct.begin(), 10000)) {
1248 19 Mar 08 peter 68      suite.err() << "test_svd: FAILED, algorithm " << algo << std::endl;
1248 19 Mar 08 peter 69      suite.err() << "reconstruction error" << std::endl;
1248 19 Mar 08 peter 70     suite.add(false);
227 01 Feb 05 jari 71   }
227 01 Feb 05 jari 72
420 02 Dec 05 jari 73   Vtranspose*=svd.V();  // Expect unity matrix
1248 19 Mar 08 peter 74   Vtranspose += 1.0;
1248 19 Mar 08 peter 75   utility::Matrix eye(Vtranspose.rows(), Vtranspose.columns(), 1.0);
1248 19 Mar 08 peter 76   for (size_t i=0; i<eye.rows() && i<eye.columns(); ++i)
1248 19 Mar 08 peter 77     eye(i,i)=2.0;
1248 19 Mar 08 peter 78   if (!suite.equal_range(eye.begin(), eye.end(), Vtranspose.begin(), 100) ) {
4200 19 Aug 22 peter 79     suite.err() << "test_svd: FAILED, algorithm " << algo
1248 19 Mar 08 peter 80                 << " V orthogonality error " << std::endl;
1248 19 Mar 08 peter 81     suite.add(false);
227 01 Feb 05 jari 82   }
227 01 Feb 05 jari 83
1121 22 Feb 08 peter 84   utility::Matrix Utranspose(svd.U());
420 02 Dec 05 jari 85   Utranspose.transpose();
420 02 Dec 05 jari 86   Utranspose*=svd.U();  // Expect unity matrix
1248 19 Mar 08 peter 87   Utranspose += 1.0;
1248 19 Mar 08 peter 88   eye = utility::Matrix(Utranspose.rows(), Utranspose.columns(), 1.0);
1248 19 Mar 08 peter 89   for (size_t i=0; i<eye.rows() && i<eye.columns(); ++i)
1248 19 Mar 08 peter 90     eye(i,i)=2.0;
1248 19 Mar 08 peter 91   if (!suite.equal_range(eye.begin(), eye.end(), Utranspose.begin(), 100) ) {
4200 19 Aug 22 peter 92     suite.err() << "test_svd: FAILED, algorithm " << algo
1248 19 Mar 08 peter 93                 << " U orthogonality error " << std::endl;
1248 19 Mar 08 peter 94     suite.add(false);
227 01 Feb 05 jari 95   }
227 01 Feb 05 jari 96 }
227 01 Feb 05 jari 97
227 01 Feb 05 jari 98
4367 12 Sep 23 peter 99 void is_eye(const utility::Matrix& X, test::Suite& suite, bool expect)
4367 12 Sep 23 peter 100 {
4367 12 Sep 23 peter 101   if (X.rows() != X.columns()) {
4367 12 Sep 23 peter 102     suite.err() << "not square input\n";
4367 12 Sep 23 peter 103     suite.add(false);
4367 12 Sep 23 peter 104     return;
4367 12 Sep 23 peter 105   }
227 01 Feb 05 jari 106
4367 12 Sep 23 peter 107   const double margin = 1e-10;
4367 12 Sep 23 peter 108   size_t n = X.rows();
4367 12 Sep 23 peter 109   utility::Matrix I(n, n);
4367 12 Sep 23 peter 110   for (size_t i=0; i<n; ++i)
4367 12 Sep 23 peter 111     I(i,i) = 1.0;
4367 12 Sep 23 peter 112
4367 12 Sep 23 peter 113   if (X.equal(I, margin) == expect)
4367 12 Sep 23 peter 114     return;
4367 12 Sep 23 peter 115
4367 12 Sep 23 peter 116   suite.add(false);
4367 12 Sep 23 peter 117   suite.err() << "failed\n";
4367 12 Sep 23 peter 118
4367 12 Sep 23 peter 119   if (expect) {
4367 12 Sep 23 peter 120     for (size_t i=0; i<n; ++i)
4367 12 Sep 23 peter 121       for (size_t j=0; j<n; ++j)
4367 12 Sep 23 peter 122         if (!suite.equal_fix(X(i, j), I(i,j), margin))
4367 12 Sep 23 peter 123           suite.err() << "error: element: " << i << " : " << j << "\n";
4367 12 Sep 23 peter 124   }
4367 12 Sep 23 peter 125   else {
4367 12 Sep 23 peter 126     suite.err() << "expected matrix not to be I\n";
4367 12 Sep 23 peter 127   }
4367 12 Sep 23 peter 128 }
4367 12 Sep 23 peter 129
4367 12 Sep 23 peter 130
4367 12 Sep 23 peter 131 void do_test2(test::Suite& suite)
4367 12 Sep 23 peter 132 {
4367 12 Sep 23 peter 133   utility::Matrix A = test::generate_Matrix(10,4,2.5);
4367 12 Sep 23 peter 134   utility::SVD svd(A);
4367 12 Sep 23 peter 135   svd.decompose();
4367 12 Sep 23 peter 136   utility::Matrix S(svd.s().size(), svd.s().size());
4367 12 Sep 23 peter 137   for (size_t i=0; i<svd.s().size(); ++i)
4367 12 Sep 23 peter 138     S(i,i) = svd.s()(i);
4367 12 Sep 23 peter 139   utility::Matrix B = svd.U() * S * transpose(svd.V());
4367 12 Sep 23 peter 140   if (!suite.add(suite.equal_matrix(A, B, 1000)))
4367 12 Sep 23 peter 141     suite.err() << "error: U * S * V' does not equal A\n";
4367 12 Sep 23 peter 142
4367 12 Sep 23 peter 143   suite.out() << "test U * U' != I\n";
4367 12 Sep 23 peter 144   is_eye(svd.U() * transpose(svd.U()), suite, false);
4367 12 Sep 23 peter 145   suite.out() << "test U' * U == I\n";
4367 12 Sep 23 peter 146   is_eye(transpose(svd.U()) * svd.U(), suite, true);
4367 12 Sep 23 peter 147   suite.out() << "test V * V' == I\n";
4367 12 Sep 23 peter 148   is_eye(svd.V() * transpose(svd.V()), suite, true);
4367 12 Sep 23 peter 149   suite.out() << "test V' * V == I\n";
4367 12 Sep 23 peter 150   is_eye(transpose(svd.V()) * svd.V(), suite, true);
4367 12 Sep 23 peter 151 }
4367 12 Sep 23 peter 152
4367 12 Sep 23 peter 153
1248 19 Mar 08 peter 154 int main(int argc, char* argv[])
227 01 Feb 05 jari 155 {
1248 19 Mar 08 peter 156   theplu::yat::test::Suite suite(argc, argv);
227 01 Feb 05 jari 157
227 01 Feb 05 jari 158   // The GSL Jacobi, Golub-Reinsch, and modified Golub-Reinsch
227 01 Feb 05 jari 159   // implementations supports rows>=columns matrix dimensions only
1248 19 Mar 08 peter 160   do_test(12,12,utility::SVD::GolubReinsch, suite);
1248 19 Mar 08 peter 161   do_test(12,4,utility::SVD::GolubReinsch, suite);
1248 19 Mar 08 peter 162   do_test(12,12,utility::SVD::ModifiedGolubReinsch, suite);
1248 19 Mar 08 peter 163   do_test(12,4,utility::SVD::ModifiedGolubReinsch, suite);
1248 19 Mar 08 peter 164   do_test(12,12,utility::SVD::Jacobi, suite);
1248 19 Mar 08 peter 165   do_test(12,4,utility::SVD::Jacobi, suite);
4367 12 Sep 23 peter 166   do_test2(suite);
227 01 Feb 05 jari 167
1248 19 Mar 08 peter 168   return suite.return_value();
227 01 Feb 05 jari 169 }

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