Gauss-Seidel iterative method for Kernel Ridge Regression learning

shogun-toolbox · Apr 18, 2012 · 0f16159 · 0f16159
1 parent 024fc0a
commit 0f16159
Show file tree

Hide file tree

Showing 2 changed files with 144 additions and 4 deletions.
diff --git a/src/shogun/regression/KernelRidgeRegression.cpp b/src/shogun/regression/KernelRidgeRegression.cpp
@@ -24,14 +24,16 @@ CKernelRidgeRegression::CKernelRidgeRegression()
 	init();
 }
 
-CKernelRidgeRegression::CKernelRidgeRegression(float64_t tau, CKernel* k, CLabels* lab)
+CKernelRidgeRegression::CKernelRidgeRegression(float64_t tau, CKernel* k, CLabels* lab, ETrainingType m)
 : CKernelMachine()
 {
 	init();
 
 	m_tau=tau;
 	set_labels(lab);
 	set_kernel(k);
+    set_epsilon(0.0001);
+    m_train_func=getTrainFunction(m);
 }
 
 void CKernelRidgeRegression::init()
@@ -41,7 +43,23 @@ void CKernelRidgeRegression::init()
 	SG_ADD(&m_tau, "tau", "Regularization parameter", MS_AVAILABLE);
 }
 
-bool CKernelRidgeRegression::train_machine(CFeatures* data)
+CKernelRidgeRegression::train_method CKernelRidgeRegression::getTrainFunction(ETrainingType id)
+{
+    switch (id)
+    {
+        case PINV:
+            return &CKernelRidgeRegression::train_machine_pinv;
+        break;
+        case GS:
+            return &CKernelRidgeRegression::train_machine_gs;
+        break;
+        default:
+            return &CKernelRidgeRegression::train_machine_pinv;
+        break;
+    }
+}
+
+bool CKernelRidgeRegression::train_machine_pinv(CFeatures* data)
 {
 	if (!m_labels)
 		SG_ERROR("No labels set\n");
@@ -86,6 +104,80 @@ bool CKernelRidgeRegression::train_machine(CFeatures* data)
 	return true;
 }
 
+bool CKernelRidgeRegression::train_machine_gs(CFeatures* data)
+{
+    if (!m_labels)
+        SG_ERROR("No labels set\n");
+
+    if (data)
+    {
+        if (m_labels->get_num_labels() != data->get_num_vectors())
+            SG_ERROR("Number of training vectors does not match number of labels\n");
+        kernel->init(data, data);
+    }
+    ASSERT(kernel && kernel->has_features());
+
+    // Get kernel matrix
+    SGMatrix<float64_t> A = kernel->get_kernel_matrix<float64_t>();
+    int32_t n = A.num_cols;
+    int32_t m = A.num_rows;
+    ASSERT(A.matrix && m>0 && n>0);
+
+    for(int32_t i=0; i < n; i++)
+        A.matrix[i+i*n]+=m_tau;
+
+    // re-set alphas of kernel machine
+    m_alpha.destroy_vector();
+    SGVector<float64_t> b;
+    float64_t alpha_old;
+
+    b=m_labels->get_labels_copy();
+    m_alpha=m_labels->get_labels_copy();
+    m_alpha.zero();
+
+    // tell kernel machine that all alphas are needed as 'support vectors'
+    m_svs.destroy_vector();
+    m_svs=SGVector<index_t>(m_alpha.vlen);
+    m_svs.range_fill();
+
+    if (get_alphas().vlen!=n)
+    {
+        SG_ERROR("Number of labels does not match number of kernel"
+                " columns (num_labels=%d cols=%d\n", m_alpha.vlen, n);
+    }
+
+    // Gauss-Seidel iterative method
+    float64_t sigma, err, d;
+    bool flag=true;
+    while(flag)
+    {
+        err=0.0;
+        for(int32_t i=0; i<n; i++)
+        {
+            sigma=b[i];
+            for(int32_t j=0; j<n; j++)
+                if (i!=j)
+                    sigma-=A.matrix[j+i*n]*m_alpha[j];
+            alpha_old=m_alpha[i];
+            m_alpha[i]=sigma/A.matrix[i+i*n];
+            d=fabs(alpha_old-m_alpha[i]);
+            if(d>err)
+                err=d;
+        }
+        if (err<=m_epsilon)
+            flag=false;
+    }
+
+    SG_FREE(A.matrix);
+
+    return true;
+}
+
+bool CKernelRidgeRegression::train_machine(CFeatures *data)
+{
+    return (this->*m_train_func)(data);
+}
+
 bool CKernelRidgeRegression::load(FILE* srcfile)
 {
 	SG_SET_LOCALE_C;

diff --git a/src/shogun/regression/KernelRidgeRegression.h b/src/shogun/regression/KernelRidgeRegression.h
@@ -57,21 +57,40 @@ class CKernelRidgeRegression : public CKernelMachine
 		/** default constructor */
 		CKernelRidgeRegression();
 
+        typedef bool (CKernelRidgeRegression::*train_method)(CFeatures* data);
+
+        enum ETrainingType
+        {
+            PINV=1,
+            GS=2,
+        };
+
 		/** constructor
 		 *
 		 * @param tau regularization constant tau
 		 * @param k kernel
 		 * @param lab labels
 		 */
-		CKernelRidgeRegression(float64_t tau, CKernel* k, CLabels* lab);
+        CKernelRidgeRegression(float64_t tau, CKernel* k, CLabels* lab, ETrainingType m=PINV);
 		virtual ~CKernelRidgeRegression() {}
 
+        /** set training function
+         *
+         */
+        inline void set_train_func(train_method m) { m_train_func = m; };
+
 		/** set regularization constant
 		 *
 		 * @param tau new tau
 		 */
 		inline void set_tau(float64_t tau) { m_tau = tau; };
 
+        /** set precision
+         *
+         * @param tau new tau
+         */
+        inline void set_epsilon(float64_t epsilon) { m_epsilon = epsilon; }
+
 		/** load regression from file
 		 *
 		 * @param srcfile file to load from
@@ -96,7 +115,7 @@ class CKernelRidgeRegression : public CKernelMachine
 		}
 
 		/** @return object name */
-		inline virtual const char* get_name() const { return "KernelRidgeRegression"; }
+        inline virtual const char* get_name() const { return "KernelRidgeRegression"; }
 
 	protected:
 		/** train regression
@@ -109,13 +128,42 @@ class CKernelRidgeRegression : public CKernelMachine
 		 */
 		virtual bool train_machine(CFeatures* data=NULL);
 
+        /** train regression using Gauss-Seidel iterative method
+          *
+          * @param data training data (parameter can be avoided if distance or
+          * kernel-based regressors are used and distance/kernels are
+          * initialized with train data)
+          *
+          * @return whether training was successful
+          */
+        bool train_machine_gs(CFeatures* data=NULL);
+
+        /** train regression using pinv
+          *
+          * @param data training data (parameter can be avoided if distance or
+          * kernel-based regressors are used and distance/kernels are
+          * initialized with train data)
+          *
+          * @return whether training was successful
+          */
+        bool train_machine_pinv(CFeatures* data=NULL);
+
 	private:
 		void init();
 
 	private:
 		/** regularization parameter tau */
 		float64_t m_tau;
+
+        /** epsilon constant */
+        float64_t m_epsilon;
+
+        /** training function */
+        train_method m_train_func;
+
+        static train_method getTrainFunction(ETrainingType id);
 };
 }
+
 #endif // HAVE_LAPACK
 #endif // _KERNELRIDGEREGRESSION_H__