yat/utility/KernelPCA.cc

Code
Comments
Other
Rev Date Author Line
2324 21 Sep 10 peter 1 // $Id$
2324 21 Sep 10 peter 2
2324 21 Sep 10 peter 3 /*
4359 23 Aug 23 peter 4   Copyright (C) 2010, 2012, 2022, 2023 Peter Johansson
2324 21 Sep 10 peter 5
2324 21 Sep 10 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
2324 21 Sep 10 peter 7
2324 21 Sep 10 peter 8   The yat library is free software; you can redistribute it and/or
2324 21 Sep 10 peter 9   modify it under the terms of the GNU General Public License as
2324 21 Sep 10 peter 10   published by the Free Software Foundation; either version 3 of the
2324 21 Sep 10 peter 11   License, or (at your option) any later version.
2324 21 Sep 10 peter 12
2324 21 Sep 10 peter 13   The yat library is distributed in the hope that it will be useful,
2324 21 Sep 10 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
2324 21 Sep 10 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2324 21 Sep 10 peter 16   General Public License for more details.
2324 21 Sep 10 peter 17
2324 21 Sep 10 peter 18   You should have received a copy of the GNU General Public License
2324 21 Sep 10 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
2324 21 Sep 10 peter 20 */
2324 21 Sep 10 peter 21
2881 18 Nov 12 peter 22 #include <config.h>
2881 18 Nov 12 peter 23
2324 21 Sep 10 peter 24 #include "KernelPCA.h"
2324 21 Sep 10 peter 25
2324 21 Sep 10 peter 26 #include "Matrix.h"
2324 21 Sep 10 peter 27 #include "PCA.h"
2324 21 Sep 10 peter 28
2324 21 Sep 10 peter 29 #include <gsl/gsl_linalg.h>
2324 21 Sep 10 peter 30
2324 21 Sep 10 peter 31 #include <cassert>
2324 21 Sep 10 peter 32 #include <cmath>
2324 21 Sep 10 peter 33
2324 21 Sep 10 peter 34 namespace theplu {
2324 21 Sep 10 peter 35 namespace yat {
2324 21 Sep 10 peter 36 namespace utility {
2324 21 Sep 10 peter 37
2324 21 Sep 10 peter 38   class KernelPCA::Impl
2324 21 Sep 10 peter 39   {
2324 21 Sep 10 peter 40   public:
4125 14 Jan 22 peter 41     Impl(const MatrixBase& kernel)
2324 21 Sep 10 peter 42       : projection_(0,0)
2324 21 Sep 10 peter 43     {
2324 21 Sep 10 peter 44       assert(kernel.rows()==kernel.columns());
2324 21 Sep 10 peter 45
2324 21 Sep 10 peter 46       // decompose K = Z' * Z
2324 21 Sep 10 peter 47
2324 21 Sep 10 peter 48       // k is typically based on centralized data, which implies k has
2324 21 Sep 10 peter 49       // one zero eigenvalue, i.e., the samples span a hyperspace
2324 21 Sep 10 peter 50       // N-1. This might cause numerical problems so therefore
2324 21 Sep 10 peter 51       // translate the data points away (perpendicular) from the
2324 21 Sep 10 peter 52       // hyperspace. In practice that means that we add all
2324 21 Sep 10 peter 53       // scalar-products with 1.0.
2324 21 Sep 10 peter 54       //      k += 1.0;
2324 21 Sep 10 peter 55
2324 21 Sep 10 peter 56       Matrix Z(kernel);
2324 21 Sep 10 peter 57       Z += 1.0e-10;
2324 21 Sep 10 peter 58       gsl_linalg_cholesky_decomp(Z.gsl_matrix_p());
4352 21 Aug 23 peter 59       // Set triangle below diagonal to zero
2324 21 Sep 10 peter 60       for (size_t row=0; row<Z.rows(); ++row)
2324 21 Sep 10 peter 61         for (size_t column=0; column<row; ++column)
2324 21 Sep 10 peter 62           Z(row, column) = 0.0;
2324 21 Sep 10 peter 63
2324 21 Sep 10 peter 64       PCA pca(Z);
3519 05 Oct 16 peter 65
2324 21 Sep 10 peter 66       eigenvalues_ = pca.eigenvalues();
2324 21 Sep 10 peter 67       projection_ = pca.projection(Z);
2324 21 Sep 10 peter 68     }
2324 21 Sep 10 peter 69
2324 21 Sep 10 peter 70     Vector eigenvalues_;
2324 21 Sep 10 peter 71     Matrix projection_;
2324 21 Sep 10 peter 72   };
2324 21 Sep 10 peter 73
2324 21 Sep 10 peter 74
4125 14 Jan 22 peter 75   KernelPCA::KernelPCA(const MatrixBase& kernel)
2324 21 Sep 10 peter 76     : impl_(new KernelPCA::Impl(kernel))
2324 21 Sep 10 peter 77   {
2324 21 Sep 10 peter 78   }
2324 21 Sep 10 peter 79
2324 21 Sep 10 peter 80
2324 21 Sep 10 peter 81   KernelPCA::~KernelPCA(void)
2324 21 Sep 10 peter 82   {
2324 21 Sep 10 peter 83     delete impl_;
2324 21 Sep 10 peter 84     impl_ = NULL;
2324 21 Sep 10 peter 85   }
2324 21 Sep 10 peter 86
2324 21 Sep 10 peter 87
2324 21 Sep 10 peter 88   const Vector& KernelPCA::eigenvalues(void) const
2324 21 Sep 10 peter 89   {
2324 21 Sep 10 peter 90     return impl_->eigenvalues_;
2324 21 Sep 10 peter 91   }
2324 21 Sep 10 peter 92
2324 21 Sep 10 peter 93
2324 21 Sep 10 peter 94   const Matrix& KernelPCA::projection(void) const
2324 21 Sep 10 peter 95   {
2324 21 Sep 10 peter 96     return impl_->projection_;
2324 21 Sep 10 peter 97   }
2324 21 Sep 10 peter 98
2324 21 Sep 10 peter 99 }}} // of namespace utility, yat, and theplu

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