Merge pull request #3172 from cfjhallgren/nystrom_krr

Implement the Nystrom approximate algorithm for KRR

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 <limits>
+#include <shogun/regression/KRRNystrom.h>
+#include <shogun/labels/RegressionLabels.h>
+#include <shogun/mathematics/eigen3.h>
+#include <shogun/mathematics/Math.h>
+
+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<int32_t> CKRRNystrom::subsample_indices()
+{
+	int32_t n=kernel->get_num_vec_lhs();
+	SGVector<int32_t> temp(n);
+	temp.range_fill();
+	CMath::permute(temp);
+	SGVector<int32_t> col(m_num_rkhs_basis);
+	for (auto i=0; i<m_num_rkhs_basis; ++i)
+		col[i]=temp[i];
+	CMath::qsort(col.vector, m_num_rkhs_basis);
+
+	return col;
+}
+
+bool CKRRNystrom::solve_krr_system()
+{
+	int32_t n=kernel->get_num_vec_lhs();
+
+	SGVector<float64_t> y=((CRegressionLabels*)m_labels)->get_labels();
+	if (y==NULL)
+		SG_ERROR("Labels not set.\n");
+	SGVector<int32_t> col=subsample_indices();
+	SGMatrix<float64_t> K_mm(m_num_rkhs_basis, m_num_rkhs_basis);
+	SGMatrix<float64_t> K_nm(n, m_num_rkhs_basis);
+	#pragma omp parallel for
+	for (index_t j=0; j<m_num_rkhs_basis; ++j)
+	{
+		for (index_t i=0; i<n; ++i)
+			K_nm(i,j)=kernel->kernel(i,col[j]);
+	}
+	#pragma omp parallel for
+	for (index_t i=0; i<m_num_rkhs_basis; ++i)
+		memcpy(K_mm.matrix+i*m_num_rkhs_basis, K_nm.get_row_vector(col[i]), m_num_rkhs_basis*sizeof(float64_t));
+
+	Map<MatrixXd> K_mm_eig(K_mm.matrix, m_num_rkhs_basis, m_num_rkhs_basis);
+	Map<MatrixXd> K_nm_eig(K_nm.matrix, n, m_num_rkhs_basis);
+	MatrixXd K_mn_eig = K_nm_eig.transpose();
+	Map<VectorXd> 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<MatrixXd> 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<float64_t>::epsilon();
+	const float64_t tolerance=m_num_rkhs_basis*dbl_epsilon*D.maxCoeff();
+	for (index_t i=0; i<m_num_rkhs_basis; ++i)
+	{
+		if (D(i,i)<tolerance)
+			D(i,i)=0;
+		else
+			D(i,i)=1/D(i,i);
+	}
+	MatrixXd pseudoinv=eigvec*D*eigvec.transpose();
+	alphas_eig=pseudoinv*K_mn_eig*y_eig;
+
+	/* Expand alpha with zeros to size n */
+	SGVector<float64_t> alpha_n(n);
+	alpha_n.zero();
+	for (index_t i=0; i<m_num_rkhs_basis; ++i)
+		alpha_n[col[i]]=alphas_eig[i];
+	m_alpha=alpha_n;
+
+	return true;
+}

src/shogun/regression/KRRNystrom.h

@@ -0,0 +1,128 @@
+
+/*
+ * 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.
+ */
+
+#ifndef _KRRNYSTROM_H__
+#define _KRRNYSTROM_H__
+
+#include <shogun/regression/KernelRidgeRegression.h>
+
+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<int32_t> with sampled indices
+	 */
+	SGVector<int32_t> subsample_indices();
+
+	/** Number of columns/rows to be sampled */
+	int32_t m_num_rkhs_basis;
+
+private:
+	void init();
+
+};
+
+}
+
+#endif // _KRRNYSTROM_H__

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;