Merge pull request #3044 from youssef-emad/develop

Add automatic bias computation in LinearMachine

src/shogun/classifier/svm/LibLinear.cpp

@@ -56,6 +56,8 @@ void CLibLinear::init()
 	C2=1;
 	set_max_iterations();
 	epsilon=1e-5;
+	/** Prevent default bias computation*/
+	set_compute_bias(false);
 
 	SG_ADD(&C1, "C1", "C Cost constant 1.", MS_AVAILABLE);
 	SG_ADD(&C2, "C2", "C Cost constant 2.", MS_AVAILABLE);

src/shogun/machine/LinearMachine.cpp

@@ -12,26 +12,37 @@
 #include <shogun/labels/RegressionLabels.h>
 #include <shogun/features/DotFeatures.h>
 #include <shogun/labels/Labels.h>
+#include <shogun/mathematics/eigen3.h>
 
 using namespace shogun;
+using namespace Eigen;
 
-CLinearMachine::CLinearMachine()
-: CMachine(), bias(0), features(NULL)
+CLinearMachine::CLinearMachine(): CMachine()
 {
 	init();
 }
 
-CLinearMachine::CLinearMachine(CLinearMachine* machine) : CMachine(),
-	bias(0), features(NULL)
+CLinearMachine::CLinearMachine(bool compute_bias): CMachine()
 {
-	set_w(machine->get_w().clone());
-	set_bias(machine->get_bias());
+	init();
+
+	m_compute_bias = compute_bias;
+}
 
+CLinearMachine::CLinearMachine(CLinearMachine* machine) : CMachine()
+{
 	init();
+
+	set_w(machine->get_w().clone());
+	set_bias(machine->get_bias());
 }
 
 void CLinearMachine::init()
 {
+	bias = 0;
+	features = NULL;
+	m_compute_bias = true;
+
 	SG_ADD(&w, "w", "Parameter vector w.", MS_NOT_AVAILABLE);
 	SG_ADD(&bias, "bias", "Bias b.", MS_NOT_AVAILABLE);
 	SG_ADD((CSGObject**) &features, "features", "Feature object.",
@@ -103,6 +114,16 @@ float64_t CLinearMachine::get_bias()
 	return bias;
 }
 
+void CLinearMachine::set_compute_bias(bool compute_bias)
+{
+	m_compute_bias = compute_bias;
+}
+
+bool CLinearMachine::get_compute_bias()
+{
+	return m_compute_bias;
+}
+
 void CLinearMachine::set_features(CDotFeatures* feat)
 {
 	SG_REF(feat);
@@ -120,3 +141,51 @@ void CLinearMachine::store_model_features()
 {
 }
 
+void CLinearMachine::compute_bias(CFeatures* data)
+{
+    REQUIRE(m_labels,"No labels set\n");
+
+    if (!data)
+    	data=features;
+
+    REQUIRE(data,"No features provided and no featured previously set\n");
+
+    REQUIRE(m_labels->get_num_labels() == data->get_num_vectors(),
+    	"Number of training vectors (%d) does not match number of labels (%d)\n",
+    	m_labels->get_num_labels(), data->get_num_vectors());
+
+    SGVector<float64_t> outputs = apply_get_outputs(data);
+
+    int32_t num_vec=data->get_num_vectors();
+
+    Map<VectorXd> eigen_outputs(outputs,num_vec);
+    Map<VectorXd> eigen_labels(((CRegressionLabels*)m_labels)->get_labels(),num_vec);
+
+     set_bias((eigen_labels - eigen_outputs).mean()) ;
+}
+
+
+bool CLinearMachine::train(CFeatures* data)
+{
+    /* not allowed to train on locked data */
+    if (m_data_locked)
+    {
+    	SG_ERROR("train data_lock() was called, only train_locked() is"
+    		" possible. Call data_unlock if you want to call train()\n",
+    		get_name());
+    }
+
+    if (train_require_labels())
+    {
+    	REQUIRE(m_labels,"No labels given",this->get_name());
+
+    	m_labels->ensure_valid(get_name());
+    }
+
+    bool result = train_machine(data);
+
+    if(m_compute_bias)
+    	compute_bias(data);
+
+    return result;
+}

src/shogun/machine/LinearMachine.h

@@ -66,12 +66,25 @@ class CLinearMachine : public CMachine
 		/** default constructor */
 		CLinearMachine();
 
+		/** Constructor
+		 *
+		 * @param  compute_bias new m_compute_bias
+		 * Determines if bias_compution is considered or not
+		 */
+		CLinearMachine(bool compute_bias);
+
 		/** destructor */
 		virtual ~CLinearMachine();
 
 		/** copy constructor */
 		CLinearMachine(CLinearMachine* machine);
 
+		/** Train machine
+		 *
+		 * @return whether training was successful
+		 */
+		virtual bool train(CFeatures* data=NULL);
+
 		/** get w
 		 *
 		 * @return weight vector
@@ -96,6 +109,19 @@ class CLinearMachine : public CMachine
 		 */
 		virtual float64_t get_bias();
 
+		/** Set m_compute_bias
+		 *
+		 * Determines if bias compution is considered or not
+		 * @param  compute_bias new m_compute_bias
+		 */
+		virtual void set_compute_bias(bool compute_bias);
+
+		/** Get compute bias
+		 *
+		 * @return compute_bias
+		 */
+		virtual bool get_compute_bias();		
+
 		/** set features
 		 *
 		 * @param feat features to set
@@ -149,6 +175,12 @@ class CLinearMachine : public CMachine
 		 */
 		virtual void store_model_features();
 
+		/** Computes the added bias. The bias is computed 
+		 * as the mean error between the predictions and 
+		 * the true labels.
+		*/
+		void compute_bias(CFeatures* data);
+
 	private:
 
 		void init();
@@ -160,6 +192,8 @@ class CLinearMachine : public CMachine
 		float64_t bias;
 		/** features */
 		CDotFeatures* features;
+		/** If true, bias is computed in ::train method */
+		bool m_compute_bias;
 };
 }
 #endif

src/shogun/regression/LinearRidgeRegression.cpp

@@ -6,16 +6,17 @@
  *
  * Copyright (C) 2012 Soeren Sonnenburg
  */
-
 #include <shogun/lib/config.h>
 
 #ifdef HAVE_LAPACK
 #include <shogun/regression/LinearRidgeRegression.h>
 #include <shogun/mathematics/lapack.h>
 #include <shogun/mathematics/Math.h>
 #include <shogun/labels/RegressionLabels.h>
+#include <shogun/mathematics/eigen3.h>
 
 using namespace shogun;
+using namespace Eigen;
 
 CLinearRidgeRegression::CLinearRidgeRegression()
 : CLinearMachine()
@@ -42,58 +43,60 @@ void CLinearRidgeRegression::init()
 
 bool CLinearRidgeRegression::train_machine(CFeatures* data)
 {
-	if (!m_labels)
-		SG_ERROR("No labels set\n")
-
-	if (!data)
-		data=features;
-
-	if (!data)
-		SG_ERROR("No features set\n")
+    REQUIRE(m_labels,"No labels set\n");
 
-	if (m_labels->get_num_labels() != data->get_num_vectors())
-		SG_ERROR("Number of training vectors does not match number of labels\n")
+    if (!data)
+    	data=features;
 
-	if (data->get_feature_class() != C_DENSE)
-		SG_ERROR("Expected Dense Features\n")
+    REQUIRE(data,"No features provided and no featured previously set\n");
 
-	if (data->get_feature_type() != F_DREAL)
-		SG_ERROR("Expected Real Features\n")
+    REQUIRE(m_labels->get_num_labels() == data->get_num_vectors(),
+    	"Number of training vectors (%d) does not match number of labels (%d)\n",
+    	m_labels->get_num_labels(), data->get_num_vectors());
 
-	CDenseFeatures<float64_t>* feats=(CDenseFeatures<float64_t>*) data;
-	int32_t num_feat=feats->get_num_features();
-	int32_t num_vec=feats->get_num_vectors();
+    REQUIRE(data->get_feature_class() == C_DENSE,
+    	"Expected Dense Features (%d) but got (%d)\n",
+    	C_DENSE, data->get_feature_class());
 
-	// Get kernel matrix
-	SGMatrix<float64_t> kernel_matrix(num_feat,num_feat);
-	SGVector<float64_t> y(num_feat);
+    REQUIRE(data->get_feature_type() == F_DREAL,
+    	"Expected Real Features (%d) but got (%d)\n",
+    	F_DREAL, data->get_feature_type());
 
-	// init
-	kernel_matrix.zero();
-	y.zero();
+    CDenseFeatures<float64_t>* feats=(CDenseFeatures<float64_t>*) data;
+    int32_t num_feat=feats->get_num_features();
+    int32_t num_vec=feats->get_num_vectors();
 
-	for (int32_t i=0; i<num_feat; i++)
-		kernel_matrix.matrix[i+i*num_feat]+=m_tau;
+    SGMatrix<float64_t> kernel_matrix(num_feat,num_feat);
+    SGMatrix<float64_t> feats_matrix(feats->get_feature_matrix());
+    SGVector<float64_t> y(num_feat);
+    SGVector<float64_t> tau_vector(num_feat);
 
-	for (int32_t i=0; i<num_vec; i++)
-	{
-		SGVector<float64_t> v = feats->get_feature_vector(i);
-		ASSERT(v.vlen==num_feat)
+    tau_vector.zero();
+    tau_vector.add(m_tau);
 
-		cblas_dger(CblasColMajor, num_feat,num_feat, 1.0, v.vector,1,
-				v.vector,1, kernel_matrix.matrix, num_feat);
+    Map<MatrixXd> eigen_kernel_matrix(kernel_matrix.matrix, num_feat,num_feat);
+    Map<MatrixXd> eigen_feats_matrix(feats_matrix.matrix, num_feat,num_vec);
+    Map<VectorXd> eigen_y(y.vector, num_feat);
+    Map<VectorXd> eigen_labels(((CRegressionLabels*)m_labels)->get_labels(),num_vec);
+    Map<VectorXd> eigen_tau(tau_vector.vector, num_feat);
 
-		cblas_daxpy(num_feat, ((CRegressionLabels*) m_labels)->get_label(i), v.vector, 1, y.vector, 1);
+    eigen_kernel_matrix = eigen_feats_matrix*eigen_feats_matrix.transpose();
 
-		feats->free_feature_vector(v, i);
-	}
+    eigen_kernel_matrix.diagonal() += eigen_tau;
 
-	clapack_dposv(CblasRowMajor,CblasUpper, num_feat, 1, kernel_matrix.matrix, num_feat,
-			y.vector, num_feat);
+    eigen_y = eigen_feats_matrix*eigen_labels ;
 
-	set_w(y);
+    LLT<MatrixXd> llt;
+    llt.compute(eigen_kernel_matrix);
+    if(llt.info() != Eigen::Success)
+    {
+    	SG_WARNING("Features covariance matrix was not positive definite\n");
+    	return false;
+    }
+    eigen_y = llt.solve(eigen_y);
 
-	return true;
+    set_w(y);
+    return true;
 }
 
 bool CLinearRidgeRegression::load(FILE* srcfile)

src/shogun/regression/LinearRidgeRegression.h

@@ -37,7 +37,7 @@ namespace shogun
  * {\bf w} = \left(\tau {\bf I}+ \sum_{i=1}^N{\bf x}_i{\bf x}_i^T\right)^{-1}\left(\sum_{i=1}^N y_i{\bf x}_i\right)
  * \f]
  *
- * The expressed solution is a linear method with bias 0 (cf. CLinearMachine).
+ * The expressed solution is a linear method with bias b (cf. CLinearMachine).
  */
 class CLinearRidgeRegression : public CLinearMachine
 {