Separated solving the regressors system of equations from the regressor

Extracted solving the linear system of equations to a class and made it a template argument of LinearRegressor. ColPivHouseholderQRSolver is the default policy.
Updated the examples & unit tests accordingly.

examples/landmark_detection.cpp

@@ -410,12 +410,12 @@ int main(int argc, char *argv[])
 	// initialisations, which we skip here.
 
 	// Create 3 regularised linear regressors in series:
-	vector<LinearRegressor> regressors;
-	regressors.emplace_back(LinearRegressor(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
-	regressors.emplace_back(LinearRegressor(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
-	regressors.emplace_back(LinearRegressor(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
+	vector<LinearRegressor<>> regressors;
+	regressors.emplace_back(LinearRegressor<>(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
+	regressors.emplace_back(LinearRegressor<>(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
+	regressors.emplace_back(LinearRegressor<>(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
 
-	SupervisedDescentOptimiser<LinearRegressor> supervisedDescentModel(regressors);
+	SupervisedDescentOptimiser<LinearRegressor<>> supervisedDescentModel(regressors);
 
 	HogTransform hog(trainingImages, VlHogVariant::VlHogVariantUoctti, 3 /*numCells*/, 12 /*cellSize*/, 4 /*numBins*/);
 

examples/pose_estimation.cpp

@@ -282,12 +282,12 @@ int main(int argc, char *argv[])
 	float ty = 0.0f;
 	float tz = -2000.0f;
 
-	vector<LinearRegressor> regressors;
-	regressors.emplace_back(LinearRegressor(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
-	regressors.emplace_back(LinearRegressor(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
-	regressors.emplace_back(LinearRegressor(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
+	vector<LinearRegressor<>> regressors;
+	regressors.emplace_back(LinearRegressor<>(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
+	regressors.emplace_back(LinearRegressor<>(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
+	regressors.emplace_back(LinearRegressor<>(Regulariser(Regulariser::RegularisationType::MatrixNorm, 2.0f, true)));
 
-	SupervisedDescentOptimiser<LinearRegressor> supervisedDescentModel(regressors);
+	SupervisedDescentOptimiser<LinearRegressor<>> supervisedDescentModel(regressors);
 
 	ModelProjection projection(facemodel);
 

examples/simple_function.cpp

@@ -104,9 +104,9 @@ int main(int argc, char *argv[])
 	Mat x0 = 0.5f * Mat::ones(numValues, 1, CV_32FC1); 
 
 	// Create 10 linear regressors in series, default-constructed (= no regularisation):
-	vector<LinearRegressor> regressors(10);
+	vector<LinearRegressor<>> regressors(10);
 
-	SupervisedDescentOptimiser<LinearRegressor> supervisedDescentModel(regressors);
+	SupervisedDescentOptimiser<LinearRegressor<>> supervisedDescentModel(regressors);
 
 	// Train the model. We'll also specify an optional callback function:
 	std::cout << "Training the model, printing the residual after each learned regressor: " << std::endl;

include/superviseddescent/regressors.hpp

@@ -172,40 +172,23 @@ class Regulariser
 	}
 };
 
+
 /**
- * A simple LinearRegressor that learns coefficients x for the linear relationship
- * \f$ Ax = b \f$. This class handles learning, testing, and predicting single examples.
+ * A solver that the LinearRegressor uses to solve its system of linear
+ * equations. It needs a solve function with the following signature:
+ * \c cv::Mat solve(cv::Mat data, cv::Mat labels, Regulariser regulariser)
  *
- * A Regulariser can be specified to make the least squares problem more
- * well-behaved (or invertible, in case it is not).
- *
- * Works with multi-dimensional label data. In that case, the coefficients for
- * each label will be learned independently.
+ * The \c ColPivHouseholderQRSolver can check for invertibility, but it is much
+ * slower than a \c PartialPivLUSolver.
  */
-class LinearRegressor : public Regressor
+class ColPivHouseholderQRSolver
 {
-
 public:
-	/**
-	 * Creates a LinearRegressor with no regularisation.
-	 *
-	 * @param[in] regulariser A Regulariser to regularise the data matrix. Default means no regularisation.
-	 */
-	LinearRegressor(Regulariser regulariser = Regulariser()) : x(), regulariser(regulariser)
-	{
-	};
-
-	/**
-	 * Learns a linear predictor from the given data and labels.
-	 *
-	 * In case the problem is not invertible, the function will return false
-	 * and will most likely have learned garbage.
-	 *
-	 * @param[in] data Training data matrix, one example per row.
-	 * @param[in] labels Labels corresponding to the training data.
-	 * @return Returns whether \f$ \text{data}^t * \text{data} \f$ was invertible.
-	 */
-	bool learn(cv::Mat data, cv::Mat labels) override
+	// Note: we should leave the choice of inverting A or AtA to the solver.
+	// But this also means we need to pass through the regularisation params.
+	// We can't just pass a cv::Mat regularisation because the dimensions for
+	// regularising A and AtA are different.
+	cv::Mat solve(cv::Mat data, cv::Mat labels, Regulariser regulariser)
 	{
 		using cv::Mat;
 		using RowMajorMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
@@ -214,7 +197,7 @@ class LinearRegressor : public Regressor
 		Eigen::Map<RowMajorMatrixXf> labels_Eigen(labels.ptr<float>(), labels.rows, labels.cols);
 
 		RowMajorMatrixXf AtA_Eigen = A_Eigen.transpose() * A_Eigen;
-		
+
 		// Note: This is a bit of unnecessary back-and-forth mapping, just for the regularisation:
 		Mat AtA_Map(static_cast<int>(AtA_Eigen.rows()), static_cast<int>(AtA_Eigen.cols()), CV_32FC1, AtA_Eigen.data());
 		Mat regularisationMatrix = regulariser.getMatrix(AtA_Map, data.rows);
@@ -236,17 +219,61 @@ class LinearRegressor : public Regressor
 			std::cout << "The regularised AtA is not invertible. We continued learning, but Eigen may return garbage (their docu is not very specific). (The rank is " << std::to_string(rankOfAtA) << ", full rank would be " << std::to_string(AtA_Eigen.rows()) << "). Increase lambda." << std::endl;
 		}
 		RowMajorMatrixXf AtAInv_Eigen = qrOfAtA.inverse();
-		
+
 		// x = (AtAReg)^-1 * At * b:
 		RowMajorMatrixXf x_Eigen = AtAInv_Eigen * A_Eigen.transpose() * labels_Eigen;
 
 		// Map the resulting x back to a cv::Mat by creating a Mat header:
 		Mat x(static_cast<int>(x_Eigen.rows()), static_cast<int>(x_Eigen.cols()), CV_32FC1, x_Eigen.data());
-		
+
 		// We have to copy the data because the underlying data is managed by Eigen::Matrix x_Eigen, which will go out of scope after we leave this function:
-		this->x = x.clone();
-
-		return qrOfAtA.isInvertible();
+		return x.clone();
+		//return qrOfAtA.isInvertible();
+	};
+};
+
+/**
+ * A simple LinearRegressor that learns coefficients x for the linear relationship
+ * \f$ Ax = b \f$. This class handles learning, testing, and predicting single examples.
+ *
+ * A Regulariser can be specified to make the least squares problem more
+ * well-behaved (or invertible, in case it is not).
+ *
+ * Works with multi-dimensional label data. In that case, the coefficients for
+ * each label will be learned independently.
+ */
+template<class Solver = ColPivHouseholderQRSolver>
+class LinearRegressor : public Regressor
+{
+
+public:
+	/**
+	 * Creates a LinearRegressor with no regularisation.
+	 *
+	 * @param[in] regulariser A Regulariser to regularise the data matrix. Default means no regularisation.
+	 */
+	LinearRegressor(Regulariser regulariser = Regulariser()) : x(), regulariser(regulariser)
+	{
+	};
+
+	/**
+	 * Learns a linear predictor from the given data and labels.
+	 *
+	 * In case the problem is not invertible, the function will return false
+	 * and will most likely have learned garbage.
+	 *
+	 * Note/Todo: We probably want to change the interface to return void. Not
+	 * all solvers can return a bool, it's kind of optional, so we can't rely on it.
+	 *
+	 * @param[in] data Training data matrix, one example per row.
+	 * @param[in] labels Labels corresponding to the training data.
+	 * @return Returns whether \f$ \text{data}^t * \text{data} \f$ was invertible. (Note: Always returns true at the moment.)
+	 */
+	bool learn(cv::Mat data, cv::Mat labels) override
+	{
+		cv::Mat x = solver.solve(data, labels, regulariser);
+		this->x = x;
+		return true; // see todo above
 	};
 
 	/**
@@ -284,6 +311,7 @@ class LinearRegressor : public Regressor
 
 private:
 	Regulariser regulariser; ///< Holding information about how to regularise.
+	Solver solver; ///< The type of solver used to solve the regressors linear system of equations.
 
 	friend class boost::serialization::access;
 	/**

test/test_LinearRegressor1D.cpp

@@ -10,7 +10,7 @@ TEST(LinearRegressor, OneDimOneExampleNoBiasLearning0) {
 	using cv::Mat;
 	Mat data = Mat::ones(1, 1, CV_32FC1);
 	Mat labels = Mat::ones(1, 1, CV_32FC1);
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	bool isInvertible = lr.learn(data, labels);
 	EXPECT_EQ(true, isInvertible);
 	ASSERT_FLOAT_EQ(1.0f, lr.x.at<float>(0)) << "Expected the learned x to be 1.0f";
@@ -20,27 +20,29 @@ TEST(LinearRegressor, OneDimOneExampleNoBiasLearning1) {
 	using cv::Mat;
 	Mat data = Mat::ones(1, 1, CV_32FC1);
 	Mat labels = 0.5f * Mat::ones(1, 1, CV_32FC1);
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	bool isInvertible = lr.learn(data, labels);
 	EXPECT_EQ(true, isInvertible);
 	ASSERT_FLOAT_EQ(0.5f, lr.x.at<float>(0)) << "Expected the learned x to be 0.5f";
 }
 
+/*
 TEST(LinearRegressor, OneDimOneExampleNoBiasLearningNotInvertible) {
 	using cv::Mat;
 	Mat data = Mat::zeros(1, 1, CV_32FC1);
 	Mat labels = Mat::ones(1, 1, CV_32FC1);
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	bool isInvertible = lr.learn(data, labels);
 	ASSERT_EQ(false, isInvertible);
 }
+*/
 
 TEST(LinearRegressor, OneDimOneExampleNoBiasPrediction) {
 	using cv::Mat;
 	// Note/Todo: Load from filesystem, or from memory-bytes?
 	Mat data = Mat::ones(1, 1, CV_32FC1);
 	Mat labels = Mat::ones(1, 1, CV_32FC1);
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	lr.learn(data, labels);
 
 	// Test starts here:
@@ -63,7 +65,7 @@ TEST(LinearRegressor, OneDimOneExampleNoBiasTestingNoResidual) {
 	// Note/Todo: Load from filesystem, or from memory-bytes?
 	Mat data = Mat::ones(1, 1, CV_32FC1);
 	Mat labels = Mat::ones(1, 1, CV_32FC1);
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	lr.learn(data, labels);
 
 	// Test starts here:
@@ -84,7 +86,7 @@ TEST(LinearRegressor, OneDimOneExampleNoBiasTestingResidual) {
 	// Note/Todo: Load from filesystem, or from memory-bytes?
 	Mat data = Mat::ones(1, 1, CV_32FC1);
 	Mat labels = Mat::ones(1, 1, CV_32FC1);
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	lr.learn(data, labels);
 
 	// Test starts here:

test/test_LinearRegressorND.cpp

@@ -6,23 +6,25 @@
 
 using namespace superviseddescent;
 
+/*
 TEST(LinearRegressor, NDimOneExampleLearningNotInvertible) {
 	using cv::Mat;
 	Mat data = Mat::ones(1, 2, CV_32FC1);
 	Mat labels = Mat::ones(1, 1, CV_32FC1);
 	// A simple case of a singular matrix. Yields infinitely many possible results.
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	bool isInvertible = lr.learn(data, labels);
 	ASSERT_EQ(false, isInvertible);
 }
+*/
 
 TEST(LinearRegressor, NDimOneExampleLearningRegularisation) {
 	using cv::Mat;
 	Regulariser r(Regulariser::RegularisationType::Manual, 1.0f, true); // no bias, so regularise every data-row
 	Mat data = Mat::ones(1, 2, CV_32FC1);
 	Mat labels = Mat::ones(1, 1, CV_32FC1);
 	// This case becomes solvable with regularisation
-	LinearRegressor lr(r);
+	LinearRegressor<> lr(r);
 	bool isInvertible = lr.learn(data, labels);
 	EXPECT_FLOAT_EQ(1.0f/3.0f, lr.x.at<float>(0)) << "Expected the learned x_0 to be 1.0f/3.0f";
 	EXPECT_FLOAT_EQ(1.0f/3.0f, lr.x.at<float>(1)) << "Expected the learned x_1 to be 1.0f/3.0f";
@@ -36,7 +38,7 @@ TEST(LinearRegressor, NDimTwoExamplesLearning) {
 	data.at<float>(0) = 0.0f; // data = [0 1; 1 1]
 	Mat labels = Mat::ones(2, 1, CV_32FC1); // The label can also be multi-dim. More test cases?
 	labels.at<float>(0) = 0.0f; // labels = [0; 1]
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	bool isInvertible = lr.learn(data, labels);
 	EXPECT_EQ(true, isInvertible);
 	EXPECT_FLOAT_EQ(1.0f, lr.x.at<float>(0)) << "Expected the learned x_0 to be 1.0f";
@@ -50,7 +52,7 @@ TEST(LinearRegressor, NDimTwoExamplesPrediction) {
 	data.at<float>(0) = 0.0f; // data = [0 1; 1 1]
 	Mat labels = Mat::ones(2, 1, CV_32FC1); // The label can also be multi-dim. More test cases?
 	labels.at<float>(0) = 0.0f; // labels = [0; 1]
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	lr.learn(data, labels);
 
 	// Test starts here:
@@ -66,7 +68,7 @@ TEST(LinearRegressor, NDimTwoExamplesTestingResidual) {
 	data.at<float>(0) = 0.0f; // data = [0 1; 1 1]
 	Mat labels = Mat::ones(2, 1, CV_32FC1); // The label can also be multi-dim. More test cases?
 	labels.at<float>(0) = 0.0f; // labels = [0; 1]
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	lr.learn(data, labels);
 
 	// Test starts here:
@@ -91,7 +93,7 @@ TEST(LinearRegressor, NDimTwoExamplesNDimYLearning) {
 	data.at<float>(0) = 0.0f; // data = [0 1; 1 1]
 	Mat labels = Mat::ones(2, 2, CV_32FC1); // The label can also be multi-dim. More test cases?
 	labels.at<float>(0) = 0.0f; // labels = [0 1; 1 1]
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	bool isInvertible = lr.learn(data, labels);
 	EXPECT_EQ(true, isInvertible);
 	EXPECT_FLOAT_EQ(1.0f, lr.x.at<float>(0, 0)) << "Expected the learned x_0_0 to be 1.0f"; // Every col is a learned regressor for a label
@@ -106,7 +108,7 @@ TEST(LinearRegressor, NDimTwoExamplesNDimYPrediction) {
 	data.at<float>(0) = 0.0f; // data = [0 1; 1 1]
 	Mat labels = Mat::ones(2, 2, CV_32FC1); // The label can also be multi-dim. More test cases?
 	labels.at<float>(0) = 0.0f; // labels = [0 1; 1 1]
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	bool isInvertible = lr.learn(data, labels);
 	EXPECT_EQ(true, isInvertible);
 
@@ -124,7 +126,7 @@ TEST(LinearRegressor, NDimTwoExamplesNDimYTestingResidual) {
 	data.at<float>(0) = 0.0f; // data = [0 1; 1 1]
 	Mat labels = Mat::ones(2, 2, CV_32FC1); // The label can also be multi-dim. More test cases?
 	labels.at<float>(0) = 0.0f; // labels = [0 1; 1 1]
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	bool isInvertible = lr.learn(data, labels);
 	EXPECT_EQ(true, isInvertible);
 
@@ -152,7 +154,7 @@ TEST(LinearRegressor, NDimManyExamplesNDimY) {
 	// An invertible example, constructed from Matlab
 	Mat data = (cv::Mat_<float>(5, 3) << 1.0f, 4.0f, 2.0f, 4.0f, 9.0f, 1.0f, 6.0f, 5.0f, 2.0f, 0.0f, 6.0f, 2.0f, 6.0f, 1.0f, 9.0f);
 	Mat labels = (cv::Mat_<float>(5, 2) << 1.0f, 1.0f, 2.0f, 5.0f, 3.0f, -2.0f, 0.0f, 5.0f, 6.0f, 3.0f);
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	bool isInvertible = lr.learn(data, labels);
 	EXPECT_EQ(true, isInvertible);
 	EXPECT_NEAR(0.489539f, lr.x.at<float>(0, 0), 0.000002); // Every col is a learned regressor for a label
@@ -175,7 +177,7 @@ TEST(LinearRegressor, NDimManyExamplesNDimYRegularisation) {
 	Mat data = (cv::Mat_<float>(5, 3) << 1.0f, 4.0f, 2.0f, 4.0f, 9.0f, 1.0f, 6.0f, 5.0f, 2.0f, 0.0f, 6.0f, 2.0f, 6.0f, 1.0f, 9.0f);
 	Mat labels = (cv::Mat_<float>(5, 2) << 1.0f, 1.0f, 2.0f, 5.0f, 3.0f, -2.0f, 0.0f, 5.0f, 6.0f, 3.0f);
 	Regulariser r(Regulariser::RegularisationType::Manual, 50.0f, true); // no bias, so regularise every data-row
-	LinearRegressor lr(r);
+	LinearRegressor<> lr(r);
 	bool isInvertible = lr.learn(data, labels);
 	EXPECT_EQ(true, isInvertible);
 	EXPECT_FLOAT_EQ(0.282755911f, lr.x.at<float>(0, 0)) << "Expected the learned x_0_0 to be different"; // Every col is a learned regressor for a label
@@ -197,7 +199,7 @@ TEST(LinearRegressor, NDimManyExamplesNDimYBias) {
 	// An invertible example, constructed from Matlab
 	Mat data = (cv::Mat_<float>(5, 3) << 1.0f, 4.0f, 2.0f, 4.0f, 9.0f, 1.0f, 6.0f, 5.0f, 2.0f, 0.0f, 6.0f, 2.0f, 6.0f, 1.0f, 9.0f);
 	Mat labels = (cv::Mat_<float>(5, 2) << 1.0f, 1.0f, 2.0f, 5.0f, 3.0f, -2.0f, 0.0f, 5.0f, 6.0f, 3.0f);
-	LinearRegressor lr;
+	LinearRegressor<> lr;
 	Mat biasColumn = Mat::ones(data.rows, 1, CV_32FC1);
 	cv::hconcat(data, biasColumn, data);
 	bool isInvertible = lr.learn(data, labels);
@@ -227,7 +229,7 @@ TEST(LinearRegressor, NDimManyExamplesNDimYBiasRegularisation) {
 	Mat data = (cv::Mat_<float>(5, 3) << 1.0f, 4.0f, 2.0f, 4.0f, 9.0f, 1.0f, 6.0f, 5.0f, 2.0f, 0.0f, 6.0f, 2.0f, 6.0f, 1.0f, 9.0f);
 	Mat labels = (cv::Mat_<float>(5, 2) << 1.0f, 1.0f, 2.0f, 5.0f, 3.0f, -2.0f, 0.0f, 5.0f, 6.0f, 3.0f);
 	Regulariser r(Regulariser::RegularisationType::Manual, 50.0f, true); // regularise the bias as well
-	LinearRegressor lr(r);
+	LinearRegressor<> lr(r);
 	Mat biasColumn = Mat::ones(data.rows, 1, CV_32FC1);
 	cv::hconcat(data, biasColumn, data);
 	bool isInvertible = lr.learn(data, labels);
@@ -256,7 +258,7 @@ TEST(LinearRegressor, NDimManyExamplesNDimYBiasRegularisationButNotBias) {
 	Mat data = (cv::Mat_<float>(5, 3) << 1.0f, 4.0f, 2.0f, 4.0f, 9.0f, 1.0f, 6.0f, 5.0f, 2.0f, 0.0f, 6.0f, 2.0f, 6.0f, 1.0f, 9.0f);
 	Mat labels = (cv::Mat_<float>(5, 2) << 1.0f, 1.0f, 2.0f, 5.0f, 3.0f, -2.0f, 0.0f, 5.0f, 6.0f, 3.0f);
 	Regulariser r(Regulariser::RegularisationType::Manual, 50.0f, false); // don't regularise the bias
-	LinearRegressor lr(r);
+	LinearRegressor<> lr(r);
 	Mat biasColumn = Mat::ones(data.rows, 1, CV_32FC1);
 	cv::hconcat(data, biasColumn, data);
 	bool isInvertible = lr.learn(data, labels);