diff --git a/src/shogun/regression/KRRNystrom.cpp b/src/shogun/regression/KRRNystrom.cpp new file mode 100644 index 00000000000..e7d2672fca6 --- /dev/null +++ b/src/shogun/regression/KRRNystrom.cpp @@ -0,0 +1,136 @@ + +/* + * Copyright (c) The Shogun Machine Learning Toolbox + * Written (W) 2016 Fredrik Hallgren + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * The views and conclusions contained in the software and documentation are those + * of the authors and should not be interpreted as representing official policies, + * either expressed or implied, of the Shogun Development Team. + */ + +#include +#include +#include +#include +#include + +using namespace shogun; +using namespace Eigen; + +CKRRNystrom::CKRRNystrom() : CKernelRidgeRegression() +{ + init(); +} + +CKRRNystrom::CKRRNystrom(float64_t tau, int32_t m, CKernel* k, CLabels* lab) +: CKernelRidgeRegression(tau, k, lab) +{ + init(); + + m_num_rkhs_basis=m; + + int32_t n=kernel->get_num_vec_lhs(); + + REQUIRE(m_num_rkhs_basis <= n, "Number of sampled rows (%d) must be \ +less than number of data points (%d)\n", m_num_rkhs_basis, n); +} + +void CKRRNystrom::init() +{ + m_num_rkhs_basis=100; +} + +SGVector CKRRNystrom::subsample_indices() +{ + int32_t n=kernel->get_num_vec_lhs(); + SGVector temp(n); + temp.range_fill(); + CMath::permute(temp); + SGVector col(m_num_rkhs_basis); + for (auto i=0; iget_num_vec_lhs(); + + SGVector y=((CRegressionLabels*)m_labels)->get_labels(); + if (y==NULL) + SG_ERROR("Labels not set.\n"); + SGVector col=subsample_indices(); + SGMatrix K_mm(m_num_rkhs_basis, m_num_rkhs_basis); + SGMatrix K_nm(n, m_num_rkhs_basis); + #pragma omp parallel for + for (index_t j=0; jkernel(i,col[j]); + } + #pragma omp parallel for + for (index_t i=0; i K_mm_eig(K_mm.matrix, m_num_rkhs_basis, m_num_rkhs_basis); + Map K_nm_eig(K_nm.matrix, n, m_num_rkhs_basis); + MatrixXd K_mn_eig = K_nm_eig.transpose(); + Map y_eig(y.vector, n); + VectorXd alphas_eig(m_num_rkhs_basis); + + /* Calculate the Moore-Penrose pseudoinverse */ + MatrixXd Kplus=K_mn_eig*K_nm_eig+m_tau*K_mm_eig; + SelfAdjointEigenSolver solver(Kplus); + if (solver.info()!=Success) + { + SG_WARNING("Eigendecomposition failed.\n") + return false; + } + + /* Solve the system for alphas */ + MatrixXd D=solver.eigenvalues().asDiagonal(); + MatrixXd eigvec=solver.eigenvectors(); + float64_t dbl_epsilon=std::numeric_limits::epsilon(); + const float64_t tolerance=m_num_rkhs_basis*dbl_epsilon*D.maxCoeff(); + for (index_t i=0; i alpha_n(n); + alpha_n.zero(); + for (index_t i=0; i + +namespace shogun { + +/** @brief Class KRRNystrom implements the Nyström method for kernel ridge + * regression, using a low-rank approximation to the kernel matrix. + * + * The method is equivalent to ordinary kernel ridge regression, but through + * projecting the data on a subset of the data points, the full + * kernel matrix does not need to be computed, and the resulting system of + * equations for the alphas is cheaper to solve. + * + * Instead of the original linear system, the following approximate system + * is solved + * + * \f[ + * {\bf \alpha} = (\tau K_{m,m} + K_{m,n}K_{n,m})^+K_{m,n} {\bf y} + * ] + * + * where \f$K_{n,m}\f$ is a submatrix containing all n rows and the m columns + * corresponding to the m chosen training examples, \f$K_{m,n}\f$ is its + * transpose and \f$K_{m,m} is the submatrix with the m rows and columns + * corresponding to the training examples chosen. \f$+\f$ indicates the + * Moore-Penrose pseudoinverse. The complexity is \f$O(m^2n)\f$. + * + * Several ways to subsample columns/rows have been proposed. Here they are + * subsampled uniformly. To implement another sampling method one has to + * override the method 'subsample_indices'. + */ +class CKRRNystrom : public CKernelRidgeRegression +{ +public: + MACHINE_PROBLEM_TYPE(PT_REGRESSION); + + /** Default constructor */ + CKRRNystrom(); + + /** Constructor + * + * @param tau regularization parameter tau + * @param m number of rows/columns to choose + * @param k kernel + * @param lab labels + */ + CKRRNystrom(float64_t tau, int32_t m, CKernel* k, CLabels* lab); + + /** Default destructor */ + virtual ~CKRRNystrom() {} + + /** Set the number of columns/rows to choose + * + * @param m new m + */ + inline void set_num_rkhs_basis(int32_t m) + { + m_num_rkhs_basis=m; + + if (kernel!=NULL) + { + int32_t n=kernel->get_num_vec_lhs(); + + REQUIRE(m_num_rkhs_basis<=n, "Number of sampled rows (%d) must be \ +less than number of data points (%d)\n", m_num_rkhs_basis, n); + } + + }; + + /** @return object name */ + virtual const char* get_name() const { return "KRRNystrom"; } + +protected: + /** Train regression using the Nyström method. + * + * @return boolean to indicate success + */ + virtual bool solve_krr_system(); + + /** Sample indices to pick rows/columns from kernel matrix + * + * @return SGVector with sampled indices + */ + SGVector subsample_indices(); + + /** Number of columns/rows to be sampled */ + int32_t m_num_rkhs_basis; + +private: + void init(); + +}; + +} + +#endif // _KRRNYSTROM_H__ diff --git a/src/shogun/regression/KernelRidgeRegression.h b/src/shogun/regression/KernelRidgeRegression.h index 9458d8c856a..9f8995bdf8a 100644 --- a/src/shogun/regression/KernelRidgeRegression.h +++ b/src/shogun/regression/KernelRidgeRegression.h @@ -74,7 +74,7 @@ class CKernelRidgeRegression : public CKernelMachine * * @param tau new tau */ - inline void set_tau(float64_t tau) { m_tau = tau; }; + inline virtual void set_tau(float64_t tau) { m_tau = tau; }; /** set convergence precision for gauss seidel method * @@ -125,15 +125,16 @@ class CKernelRidgeRegression : public CKernelMachine * * @return boolean to indicate success */ - bool solve_krr_system(); + virtual bool solve_krr_system(); private: void init(); - private: + protected: /** regularization parameter tau */ float64_t m_tau; + private: /** epsilon constant */ float64_t m_epsilon; diff --git a/tests/unit/regression/krrnystrom_unittest.cc b/tests/unit/regression/krrnystrom_unittest.cc new file mode 100644 index 00000000000..b8770462e72 --- /dev/null +++ b/tests/unit/regression/krrnystrom_unittest.cc @@ -0,0 +1,159 @@ +/* +* Copyright (c) The Shogun Machine Learning Toolbox +* Written (w) 2016 Fredrik Hallgren +* All rights reserved. +* +* Redistribution and use in source and binary forms, with or without +* modification, are permitted provided that the following conditions are met: +* +* 1. Redistributions of source code must retain the above copyright notice, this +* list of conditions and the following disclaimer. +* 2. Redistributions in binary form must reproduce the above copyright notice, +* this list of conditions and the following disclaimer in the documentation +* and/or other materials provided with the distribution. +* +* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR +* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +* +* The views and conclusions contained in the software and documentation are those +* of the authors and should not be interpreted as representing official policies, +* either expressed or implied, of the Shogun Development Team. +*/ + +#include +#include +#include +#include +#include +#include +#include + +#ifdef HAVE_CXX11 + +using namespace shogun; + +/** + * Test the main algorithm by comparison of the alphas and the predictions, + * using all columns in the approximation to make sure results are very similar + */ +TEST(KRRNystrom, apply_and_compare_to_KRR_with_all_columns) +{ + /* data matrix dimensions */ + index_t num_vectors=10; + index_t num_features=1; + + /* training label data */ + SGVector lab(num_vectors); + + /* fill data matrix and labels */ + SGMatrix train_dat(num_features, num_vectors); + SGMatrix test_dat(num_features, num_vectors); + for (index_t i=0; i>(train_dat); + auto features_krr=some>(train_dat); + + /* testing features */ + auto test_features=some>(test_dat); + + /* training labels */ + auto labels=some(lab); + auto labels_krr=some(lab); + + /* kernel */ + auto kernel=some(features, features, 10, 0.5); + auto kernel_krr=some(features_krr, features_krr, 10, 0.5); + + /* kernel ridge regression and the nystrom approximation */ + float64_t tau=0.01; + auto nystrom=some(tau, num_vectors, kernel, labels); + auto krr=some(tau, kernel_krr, labels_krr); + + nystrom->train(); + krr->train(); + + SGVector alphas=nystrom->get_alphas(); + SGVector alphas_krr=krr->get_alphas(); + + for (index_t i=0; iapply_regression(test_features); + auto result_krr=krr->apply_regression(test_features); + + for (index_t i=0; iget_label(i), result_krr->get_label(i), 1E-5); +} + +/** + * Test the main algorithm by comparison of the alphas and the predictions, + * using a subset of the columns. + */ +TEST(KRRNystrom, apply_and_compare_to_KRR_with_column_subset) +{ + /* data matrix dimensions */ + index_t num_vectors=100; + index_t num_features=1; + index_t num_basis_rkhs=50; + + /* training label data */ + SGVector lab(num_vectors); + + /* fill data matrix and labels */ + SGMatrix train_dat(num_features, num_vectors); + SGMatrix test_dat(num_features, num_vectors); + for (index_t i=0; i>(train_dat); + auto features_krr=some>(train_dat); + + /* testing features */ + auto test_features=some>(test_dat); + + /* training labels */ + auto labels=some(lab); + auto labels_krr=some(lab); + + /* kernel */ + auto kernel=some(features, features, 10, 0.5); + auto kernel_krr=some(features_krr, features_krr, 10, 0.5); + + /* kernel ridge regression and the nystrom approximation */ + float64_t tau=0.01; + auto nystrom=some(tau, num_basis_rkhs, kernel, labels); + auto krr=some(tau, kernel_krr, labels_krr); + + nystrom->train(); + krr->train(); + + auto result=nystrom->apply_regression(test_features); + auto result_krr=krr->apply_regression(test_features); + + for (index_t i=0; iget_label(i), result_krr->get_label(i), 1E-1); +} + +#endif /* HAVE_CXX11 */