Merge branch 'master' of https://github.com/mlpack/mlpack into embedding

HISTORY.md

@@ -11,6 +11,11 @@
   * Additional functionality for the ARFF loader (#2486); use case sensitive
     categories (#2516).
 
+  * Add `bayesian_linear_regression` binding for the command-line, Python,
+    Julia, and Go.  Also called "Bayesian Ridge", this is equivalent to a
+    version of linear regression where the regularization parameter is
+    automatically tuned (#2030).
+
 ### mlpack 3.3.2
 ###### 2020-06-18
   * Added Noisy DQN to q_networks (#2446).

src/mlpack/methods/CMakeLists.txt

@@ -6,6 +6,7 @@ set(DIRS
   ann
   approx_kfn
   bias_svd
+  bayesian_linear_regression
   block_krylov_svd
   cf
   dbscan

src/mlpack/methods/bayesian_linear_regression/CMakeLists.txt

@@ -0,0 +1,21 @@
+# Define the files we need to compile
+# Anything not in this list will not be compiled into the output library
+set(SOURCES
+  bayesian_linear_regression.hpp
+  bayesian_linear_regression_impl.hpp
+  bayesian_linear_regression.cpp
+)
+
+# add directory name to sources
+set(DIR_SRCS)
+foreach(file ${SOURCES})
+  set(DIR_SRCS ${DIR_SRCS} ${CMAKE_CURRENT_SOURCE_DIR}/${file})
+endforeach()
+# append sources (with directory name) to list of all mlpack sources (used at the parent scope)
+set(MLPACK_SRCS ${MLPACK_SRCS} ${DIR_SRCS} PARENT_SCOPE)
+
+add_cli_executable(bayesian_linear_regression)
+add_python_binding(bayesian_linear_regression)
+add_julia_binding(bayesian_linear_regression)
+add_go_binding(bayesian_linear_regression)
+add_markdown_docs(bayesian_linear_regression "cli;python;julia;go" "regression")

src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression.cpp

@@ -0,0 +1,188 @@
+/**
+ * @file methods/bayesian_linear_regression/bayesian_linear_regression.cpp
+ * @author Clement Mercier 
+ *
+ * Implementation of Bayesian linear regression.
+ *
+ * mlpack is free software; you may redistribute it and/or modify it under the
+ * terms of the 3-clause BSD license.  You should have received a copy of the
+ * 3-clause BSD license along with mlpack.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#include "bayesian_linear_regression.hpp"
+#include <mlpack/core/util/log.hpp>
+#include <mlpack/core/util/timers.hpp>
+
+using namespace mlpack;
+using namespace mlpack::regression;
+
+BayesianLinearRegression::BayesianLinearRegression(const bool centerData,
+                                                   const bool scaleData,
+                                                   const size_t nIterMax,
+                                                   const double tol) :
+  centerData(centerData),
+  scaleData(scaleData),
+  nIterMax(nIterMax),
+  tol(tol),
+  responsesOffset(0.0),
+  alpha(0.0),
+  beta(0.0),
+  gamma(0.0)
+{/* Nothing to do */}
+
+double BayesianLinearRegression::Train(const arma::mat& data,
+                                       const arma::rowvec& responses)
+{
+  Timer::Start("bayesian_linear_regression");
+
+  arma::mat phi;
+  arma::rowvec t;
+  arma::colvec eigVal;
+  arma::mat eigVec;
+
+  // Preprocess the data. Center and scale.
+  responsesOffset = CenterScaleData(data, responses, phi, t);
+
+  if (!arma::eig_sym(eigVal, eigVec, arma::symmatu(phi * phi.t())))
+  {
+    Log::Fatal << "BayesianLinearRegression::Train(): Eigendecomposition "
+               << "of covariance failed!" << std::endl;
+  }
+
+  // Compute this quantities once and for all.
+  const arma::mat eigVecInv = inv(eigVec);
+  const arma::colvec eigVecInvPhitT = eigVecInv * phi * t.t();
+
+  // Initialize the hyperparameters and begin with an infinitely broad prior.
+  alpha = 1e-6;
+  beta =  1 / (var(t, 1) * 0.1);
+
+  unsigned short i = 0;
+  double deltaAlpha = 1.0, deltaBeta = 1.0, crit = 1.0;
+
+  while ((crit > tol) && (i < nIterMax))
+  {
+    deltaAlpha = -alpha;
+    deltaBeta = -beta;
+
+    // Update the solution.
+    omega = eigVec * diagmat(1 / (eigVal + (alpha / beta))) * eigVecInvPhitT;
+
+    // Update alpha.
+    gamma = sum(eigVal / (alpha / beta + eigVal));
+    alpha = gamma / dot(omega, omega);
+
+    // Update beta.
+    const arma::rowvec temp = t - omega.t() * phi;
+    beta = (data.n_cols - gamma) / dot(temp, temp);
+
+    // Compute the stopping criterion.
+    deltaAlpha += alpha;
+    deltaBeta += beta;
+    crit = std::abs(deltaAlpha / alpha + deltaBeta / beta);
+    i++;
+  }
+  // Compute the covariance matrix for the uncertainties later.
+  matCovariance = eigVec * diagmat(1 / (beta * eigVal + alpha)) * eigVecInv;
+
+  Timer::Stop("bayesian_linear_regression");
+
+  return RMSE(data, responses);
+}
+
+void BayesianLinearRegression::Predict(const arma::mat& points,
+                                       arma::rowvec& predictions) const
+{
+  // Center and scale the points before applying the model.
+  arma::mat matX;
+  CenterScaleDataPred(points, matX);
+  predictions = omega.t() * matX + responsesOffset;
+}
+
+void BayesianLinearRegression::Predict(const arma::mat& points,
+                                       arma::rowvec& predictions,
+                                       arma::rowvec& std) const
+{
+  // Center and scale the points before applying the model.
+  arma::mat matX;
+  CenterScaleDataPred(points, matX);
+  predictions = omega.t() * matX + responsesOffset;
+  // Compute the standard deviation for each point.
+  std = sqrt(Variance() + sum(matX % (matCovariance * matX), 0));
+}
+
+double BayesianLinearRegression::RMSE(const arma::mat& data,
+                                      const arma::rowvec& responses) const
+{
+  arma::rowvec predictions;
+  Predict(data, predictions);
+  return sqrt(mean(square(responses - predictions)));
+}
+
+double BayesianLinearRegression::CenterScaleData(const arma::mat& data,
+                                                 const arma::rowvec& responses,
+                                                 arma::mat& dataProc,
+                                                 arma::rowvec& responsesProc)
+{
+  if (!centerData && !scaleData)
+  {
+    dataProc = arma::mat(const_cast<double*>(data.memptr()), data.n_rows,
+                                             data.n_cols, false, true);
+    responsesProc = arma::rowvec(const_cast<double*>(responses.memptr()),
+                                                     responses.n_elem, false,
+                                                     true);
+  }
+
+  else if (centerData && !scaleData)
+  {
+    dataOffset = mean(data, 1);
+    responsesOffset = mean(responses);
+    dataProc = data.each_col() - dataOffset;
+    responsesProc = responses - responsesOffset;
+  }
+
+  else if (!centerData && scaleData)
+  {
+    dataScale = stddev(data, 0, 1);
+    dataProc = data.each_col() / dataScale;
+    responsesProc = arma::rowvec(const_cast<double*>(responses.memptr()),
+                                                     responses.n_elem, false,
+                                                     true);
+  }
+
+  else
+  {
+    dataOffset = mean(data, 1);
+    dataScale = stddev(data, 0, 1);
+    responsesOffset = mean(responses);
+    dataProc = (data.each_col() - dataOffset).each_col() / dataScale;
+    responsesProc = responses - responsesOffset;
+  }
+  return responsesOffset;
+}
+
+void BayesianLinearRegression::CenterScaleDataPred(
+    const arma::mat& data,
+    arma::mat& dataProc) const
+{
+  if (!centerData && !scaleData)
+  {
+    dataProc = arma::mat(const_cast<double*>(data.memptr()), data.n_rows,
+                         data.n_cols, false, true);
+  }
+
+  else if (centerData && !scaleData)
+  {
+    dataProc = data.each_col() - dataOffset;
+  }
+
+  else if (!centerData && scaleData)
+  {
+    dataProc = data.each_col() / dataScale;
+  }
+
+  else
+  {
+    dataProc = (data.each_col() - dataOffset).each_col() / dataScale;
+  }
+}

src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression.hpp

@@ -0,0 +1,285 @@
+/**
+ * @file methods/bayesian_linear_regression/bayesian_linear_regression.hpp
+ * @author Clement Mercier
+ *
+ * Definition of the BayesianRidge class, which performs the
+ * bayesian linear regression. According to the armadillo standards,
+ * all the functions consider data in column-major format.
+**/
+
+#ifndef MLPACK_METHODS_BAYESIAN_LINEAR_REGRESSION_HPP
+#define MLPACK_METHODS_BAYESIAN_LINEAR_REGRESSION_HPP
+
+#include <mlpack/prereqs.hpp>
+
+namespace mlpack {
+namespace regression {
+
+/**
+ * A Bayesian approach to the maximum likelihood estimation of the parameters 
+ * \f$ \omega \f$ of the linear regression model. The Complexity is governed by 
+ * the addition of a gaussian isotropic prior of precision \f$ \alpha \f$ over 
+ * \f$ \omega \f$:  
+ *
+ * \f[
+ * p(\omega|\alpha) = \mathcal{N}(\omega|0, \alpha^{-1}I)
+ * \f]
+ * 
+ * The optimization procedure calculates the posterior distribution of 
+ * \f$ \omega \f$ knowing the data by maximizing an approximation of the log 
+ * marginal likelihood derived from a type II maximum likelihood approximation. 
+ * The determination of \f$ alpha \f$ and of the noise precision \f$ beta \f$
+ * is part of the optimization process, leading to an automatic determination of 
+ * w. The model being entirely based on probabilty distributions, uncertainties 
+ * are available and easly computed for both the parameters and the predictions.
+ *
+ * The advantage over linear regression and ridge regression is that the 
+ * regularization is determined from all the training data alone without any
+ * require to an hold out method. 
+ *
+ * The code below is an implementation of the maximization of the evidence 
+ * function described in the section 3.5.2 of the C.Bishop book, Pattern
+ * Recognition and Machine Learning.  
+ *
+ * @code
+ * @article{MacKay91bayesianinterpolation,
+ *   author = {David J.C. MacKay},
+ *   title = {Bayesian Interpolation},
+ *   journal = {NEURAL COMPUTATION},
+ *   year = {1991},
+ *   volume = {4},
+ *   pages = {415--447}
+ * }
+ * @endcode
+ *
+ * @code
+ * @book{Bishop:2006:PRM:1162264,
+ *   author = {Bishop, Christopher M.},
+ *   title = {Pattern Recognition and Machine Learning (Information Science 
+ *            and Statistics)},
+ *   chapter = {3}
+ *   year = {2006},
+ *   isbn = {0387310738},
+ *   publisher = {Springer-Verlag},
+ *   address = {Berlin, Heidelberg},
+ * } 
+ * @endcode
+ *   
+ * Example of use:
+ *
+ * @code
+ * arma::mat xTrain; // Train data matrix. Column-major.
+ * arma::rowvec yTrain; // Train target values.
+ 
+ * // Train the model. Regularization strength is optimally tunned with the
+ * // training data alone by applying the Train method.
+ * BayesianLinearRegression estimator(); // Instanciate the estimator with default option.
+ * estimator.Train(xTrain, yTrain);
+ 
+ * // Prediction on test points.
+ * arma::mat xTest; // Test data matrix. Column-major.
+ * arma::rowvec predictions;
+ 
+ * estimator.Predict(xTest, prediction);
+ 
+ * arma::rowvec yTest; // Test target values.
+ * estimator.RMSE(xTest, yTest); // Evaluate using the RMSE score.
+ 
+ * // Compute the standard deviations of the predictions.
+ * arma::rowvec stds;
+ * estimator.Predict(xTest, responses, stds)
+ * @endcode
+ */
+class BayesianLinearRegression
+{
+ public:
+  /**
+   * Set the parameters of Bayesian Ridge regression object. The
+   * regularization parameter is automatically set to its optimal value by
+   * maximization of the marginal likelihood.
+   *
+   * @param centerData Whether or not center the data according to the
+   *    examples.
+   * @param scaleData Whether or not scale the data according to the
+   *    standard deviation of each feature.
+   * @param nIterMax Maximum number of iterations for convergency.
+   * @param tol Level from which the solution is considered sufficientlly 
+   *    stable.  
+   */
+  BayesianLinearRegression(const bool centerData = true,
+                           const bool scaleData = false,
+                           const size_t nIterMax = 50,
+                           const double tol = 1e-4);
+
+  /**
+   * Run BayesianLinearRegression. The input matrix (like all mlpack matrices) should be
+   * column-major -- each column is an observation and each row is a dimension.
+   * 
+   * @param data Column-major input data, dim(P, N).
+   * @param responses A vector of targets, dim(N).
+   * @return Root mean squared error.
+   */
+  double Train(const arma::mat& data,
+               const arma::rowvec& responses);
+
+  /**
+   * Predict \f$y_{i}\f$ for each data point in the given data matrix using the
+   * currently-trained Bayesian Ridge model.
+   *
+   * @param points The data points to apply the model.
+   * @param predictions y, Contains the  predicted values on completion.
+   * @return Root mean squared error computed on the train set.
+   */
+  void Predict(const arma::mat& points,
+               arma::rowvec& predictions) const;
+
+  /**
+   * Predict \f$y_{i}\f$ and the standard deviation of the predictive posterior 
+   * distribution for each data point in the given data matrix, using the
+   * currently-trained Bayesian Ridge estimator.
+   *
+   * @param points The data point to apply the model.
+   * @param predictions Vector which will contain calculated values on completion.
+   * @param std Standard deviations of the predictions.
+   */
+  void Predict(const arma::mat& points,
+               arma::rowvec& predictions,
+               arma::rowvec& std) const;
+
+  /**
+   * Compute the Root Mean Square Error between the predictions returned by the
+   * model and the true responses.
+   *
+   * @param data Data points to predict
+   * @param responses A vector of targets.
+   * @return Root mean squared error.
+   **/
+  double RMSE(const arma::mat& data,
+              const arma::rowvec& responses) const;
+
+  /**
+   * Get the solution vector. 
+   *
+   * @return omega Solution vector.
+   */
+  const arma::colvec& Omega() const { return omega; }
+
+  /**
+   * Get the precision (or inverse variance) of the gaussian prior. Train()
+   * must be called before.
+   *
+   * @return \f$ \alpha \f$
+   */
+  double Alpha() const { return alpha; }
+
+  /**
+   * Get the precision (or inverse variance) beta of the model. Train() must be
+   * called before.
+   *
+   * @return \f$ \beta \f$
+   */
+  double Beta() const { return beta; }
+
+  /**
+   * Get the estimated variance. Train() must be called before.
+   *   
+   * @return 1.0 / \f$ \beta \f$
+   */
+  double Variance() const { return 1.0 / Beta(); }
+
+  /**
+   * Get the mean vector computed on the features over the training points.
+   *
+   * @return responsesOffset
+   */
+  const arma::colvec& DataOffset() const { return dataOffset; }
+
+  /**
+   * Get the vector of standard deviations computed on the features over the 
+   * training points.
+   *
+   * @return dataOffset
+   */
+  const arma::colvec& DataScale() const { return dataScale; }
+
+  /**
+   * Get the mean value of the train responses.
+   *
+   * @return responsesOffset
+   */
+  double ResponsesOffset() const { return responsesOffset; }
+
+  /**
+   * Serialize the BayesianLinearRegression model.
+   **/
+  template<typename Archive>
+  void serialize(Archive& ar, const unsigned int /* version */);
+
+ private:
+  //! Center the data if true.
+  bool centerData;
+
+  //! Scale the data by standard deviations if true.
+  bool scaleData;
+
+  //! Maximum number of iterations for convergency.
+  size_t nIterMax;
+
+  //! Level from which the solution is considered sufficientlly stable.
+  double tol;
+
+  //! Mean vector computed over the points.
+  arma::colvec dataOffset;
+
+  //! Std vector computed over the points.
+  arma::colvec dataScale;
+
+  //! Mean of the response vector computed over the points.
+  double responsesOffset;
+
+  //! Precision of the prior pdf (gaussian).
+  double alpha;
+
+  //! Noise inverse variance.
+  double beta;
+
+  //! Effective number of parameters.
+  double gamma;
+
+  //! Solution vector
+  arma::colvec omega;
+
+  //! Covariance matrix of the solution vector omega.
+  arma::mat matCovariance;
+
+  /**
+   * Center and scale the data accordind to centerData and scaleData.
+   * Allows future modifications of new points.
+   *
+   * @param data Design matrix in column-major format, dim(P, N).
+   * @param responses A vector of targets.
+   * @param dataProc Data processed, dim(P, N).
+   * @param responsesProc Responses processed, dim(N).
+   * @return reponsesOffset Mean of responses.
+   */
+  double CenterScaleData(const arma::mat& data,
+                         const arma::rowvec& responses,
+                         arma::mat& dataProc,
+                         arma::rowvec& responsesProc);
+
+  /**
+   * Center and scale the points before prediction.
+   *
+   * @param data Design matrix in column-major format, dim(P, N).
+   * @param dataProc Data processed, dim(P, N).
+   */
+  void CenterScaleDataPred(const arma::mat& data,
+                           arma::mat& dataProc) const;
+};
+} // namespace regression
+} // namespace mlpack
+
+// Include implementation of serialize.
+#include "bayesian_linear_regression_impl.hpp"
+
+#endif

src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression_impl.hpp

@@ -0,0 +1,44 @@
+/**
+ * @file methods/bayesian_linear_regression/bayesian_linear_regression_impl.hpp
+ * @author Clement Mercier
+ *
+ * Implementation of templated BayesianLinearRegression functions.
+ *
+ * mlpack is free software; you may redistribute it and/or modify it under the
+ * terms of the 3-clause BSD license.  You should have received a copy of the
+ * 3-clause BSD license along with mlpack.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef MLPACK_METHODS_BAYESIAN_LINEAR_REGRESSION_IMPL_HPP
+#define MLPACK_METHODS_BAYESIAN_LINEAR_REGRESSION_IMPL_HPP
+
+#include "bayesian_linear_regression.hpp"
+
+namespace mlpack {
+namespace regression {
+
+/**
+ * Serialize the Bayesian linear regression model.
+ */
+template<typename Archive>
+void BayesianLinearRegression::serialize(Archive& ar,
+                                         const unsigned int /* version */)
+{
+  ar & BOOST_SERIALIZATION_NVP(centerData);
+  ar & BOOST_SERIALIZATION_NVP(scaleData);
+  ar & BOOST_SERIALIZATION_NVP(nIterMax);
+  ar & BOOST_SERIALIZATION_NVP(tol);
+  ar & BOOST_SERIALIZATION_NVP(dataOffset);
+  ar & BOOST_SERIALIZATION_NVP(dataScale);
+  ar & BOOST_SERIALIZATION_NVP(responsesOffset);
+  ar & BOOST_SERIALIZATION_NVP(alpha);
+  ar & BOOST_SERIALIZATION_NVP(beta);
+  ar & BOOST_SERIALIZATION_NVP(gamma);
+  ar & BOOST_SERIALIZATION_NVP(omega);
+  ar & BOOST_SERIALIZATION_NVP(matCovariance);
+}
+
+} // namespace regression
+} // namespace mlpack
+
+#endif

src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression_main.cpp

@@ -0,0 +1,199 @@
+/**
+ * @file methods/bayesian_linear_regression/bayesian_linear_regression_main.cpp
+ * @author Clement Mercier
+ *
+ * Executable for BayesianLinearRegression.
+ *
+ * mlpack is free software; you may redistribute it and/or modify it under the
+ * terms of the 3-clause BSD license.  You should have received a copy of the
+ * 3-clause BSD license along with mlpack.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#include <mlpack/prereqs.hpp>
+#include <mlpack/core/util/io.hpp>
+#include <mlpack/core/util/mlpack_main.hpp>
+
+#include "bayesian_linear_regression.hpp"
+
+using namespace arma;
+using namespace std;
+using namespace mlpack;
+using namespace mlpack::regression;
+using namespace mlpack::util;
+
+PROGRAM_INFO("BayesianLinearRegression",
+    // Short description.
+    "An implementation of the bayesian linear regression.",
+    // Long description.
+    "An implementation of the bayesian linear regression."
+    "\n"
+    "This model is a probabilistic view and implementation of the linear "
+    "regression. The final solution is obtained by computing a posterior "
+    "distribution from gaussian likelihood and a zero mean gaussian isotropic "
+    " prior distribution on the solution. "
+    "\n"
+    "Optimization is AUTOMATIC and does not require cross validation. "
+    "The optimization is performed by maximization of the evidence function. "
+    "Parameters are tuned during the maximization of the marginal likelihood. "
+    "This procedure includes the Ockham's razor that penalizes over complex "
+    "solutions. "
+    "\n\n"
+    "This program is able to train a Bayesian linear regression model or load "
+    "a model from file, output regression predictions for a test set, and save "
+    "the trained model to a file."
+    "\n\n"
+    "To train a BayesianLinearRegression model, the " +
+    PRINT_PARAM_STRING("input") + " and " + PRINT_PARAM_STRING("responses") +
+    "parameters must be given. The " + PRINT_PARAM_STRING("center") +
+    "and " + PRINT_PARAM_STRING("scale") + " parameters control the "
+    "centering and the normalizing options. A trained model can be saved with "
+    "the " + PRINT_PARAM_STRING("output_model") + ". If no training is desired "
+    "at all, a model can be passed via the " +
+    PRINT_PARAM_STRING("input_model") + " parameter."
+    "\n\n"
+    "The program can also provide predictions for test data using either the "
+    "trained model or the given input model.  Test points can be specified "
+    "with the " + PRINT_PARAM_STRING("test") + " parameter.  Predicted "
+    "responses to the test points can be saved with the " +
+    PRINT_PARAM_STRING("predictions") + " output parameter. The "
+    "corresponding standard deviation can be save by precising the " +
+    PRINT_PARAM_STRING("stds") + " parameter."
+    "\n\n"
+    "For example, the following command trains a model on the data " +
+    PRINT_DATASET("data") + " and responses " + PRINT_DATASET("responses") +
+    "with center set to true and scale set to false (so, Bayesian "
+    "linear regression is being solved, and then the model is saved to " +
+    PRINT_MODEL("bayesian_linear_regression_model") + ":"
+    "\n\n" +
+    PRINT_CALL("bayesian_linear_regression", "input", "data", "responses",
+               "responses", "center", 1, "scale", 0, "output_model",
+               "bayesian_linear_regression_model") +
+    "\n\n"
+    "The following command uses the " +
+    PRINT_MODEL("bayesian_linear_regression_model") + " to provide predicted " +
+    " responses for the data " + PRINT_DATASET("test") + " and save those " +
+    " responses to " + PRINT_DATASET("test_predictions") + ": "
+    "\n\n" +
+    PRINT_CALL("bayesian_linear_regression", "input_model",
+               "bayesian_linear_regression_model", "test", "test",
+               "predictions", "test_predictions") +
+    "\n\n"
+    "Because the estimator computes a predictive distribution instead of "
+    "simple point estimate, the " + PRINT_PARAM_STRING("stds") + " parameter "
+    "allows to save the prediction uncertainties: "
+    "\n\n" +
+    PRINT_CALL("bayesian_linear_regression", "input_model",
+               "bayesian_linear_regression_model", "test", "test",
+               "predictions", "test_predictions", "stds", "stds"),
+    SEE_ALSO("Bayesian Interpolation",
+        "https://authors.library.caltech.edu/13792/1/MACnc92a.pdf"),
+    SEE_ALSO("Bayesian Linear Regression, Section 3.3",
+        "MLA Bishop, Christopher M. Pattern Recognition and Machine "
+        "Learning. New York :Springer, 2006, section 3.3."),
+    SEE_ALSO("mlpack::regression::BayesianLinearRegression C++ class "
+        "documentation",
+        "@doxygen/classmlpack_1_1regression_1_1BayesianLinearRegression.html"));
+
+PARAM_MATRIX_IN("input", "Matrix of covariates (X).", "i");
+
+PARAM_ROW_IN("responses", "Matrix of responses/observations (y).", "r");
+
+PARAM_MODEL_IN(BayesianLinearRegression, "input_model", "Trained "
+               "BayesianLinearRegression model to use.", "m");
+
+PARAM_MODEL_OUT(BayesianLinearRegression, "output_model", "Output "
+                "BayesianLinearRegression model.", "M");
+
+PARAM_MATRIX_IN("test", "Matrix containing points to regress on (test "
+                "points).", "t");
+
+PARAM_MATRIX_OUT("predictions", "If --test_file is specified, this "
+                  "file is where the predicted responses will be saved.", "o");
+
+PARAM_MATRIX_OUT("stds", "If specified, this is where the standard deviations "
+    "of the predictive distribution will be saved.", "u");
+
+PARAM_FLAG("center", "Center the data and fit the intercept if enabled.", "c");
+
+PARAM_FLAG("scale", "Scale each feature by their standard deviations if "
+           "enabled.", "s");
+
+static void mlpackMain()
+{
+  bool center = IO::GetParam<bool>("center");
+  bool scale = IO::GetParam<bool>("scale");
+
+  // Check parameters -- make sure everything given makes sense.
+  RequireOnlyOnePassed({"input", "input_model"}, true);
+  if (IO::HasParam("input"))
+  {
+    RequireOnlyOnePassed({"responses"}, true, "if input data is specified, "
+        "responses must also be specified");
+  }
+  ReportIgnoredParam({{"input", false }}, "responses");
+
+  RequireAtLeastOnePassed({"predictions", "output_model", "stds"}, false,
+      "no results will be saved");
+
+  // Ignore out_predictions unless test is specified.
+  ReportIgnoredParam({{"test", false}}, "predictions");
+
+  BayesianLinearRegression* bayesLinReg;
+  if (IO::HasParam("input"))
+  {
+    Log::Info << "input detected " << std::endl;
+    // Initialize the object.
+    bayesLinReg = new BayesianLinearRegression(center, scale);
+
+    // Load covariates.  We can avoid LARS transposing our data by choosing to
+    // not transpose this data (that's why we used PARAM_TMATRIX_IN).
+    mat matX = std::move(IO::GetParam<arma::mat>("input"));
+
+    // Load responses.  The responses should be a one-dimensional vector, and it
+    // seems more likely that these will be stored with one response per line
+    // (one per row). So we should not transpose upon loading.
+    arma::rowvec responses = std::move(
+        IO::GetParam<arma::rowvec>("responses"));
+
+    if (responses.n_elem != matX.n_cols)
+    {
+      delete bayesLinReg;
+      Log::Fatal << "Number of responses must be equal to number of rows of X!"
+                 << endl;
+    }
+
+    arma::rowvec predictionsTrain;
+    // The Train method is ready to take data in column-major format.
+    bayesLinReg->Train(matX, responses);
+  }
+  else // We must have --input_model_file.
+  {
+    bayesLinReg = IO::GetParam<BayesianLinearRegression*>("input_model");
+  }
+
+  if (IO::HasParam("test"))
+  {
+    Log::Info << "Regressing on test points." << endl;
+    // Load test points.
+    mat testPoints = std::move(IO::GetParam<arma::mat>("test"));
+    arma::rowvec predictions;
+
+    if (IO::HasParam("stds"))
+    {
+      arma::rowvec std;
+      bayesLinReg->Predict(testPoints, predictions, std);
+
+      // Save the standard deviation of the test points (one per line).
+      IO::GetParam<arma::mat>("stds") = std::move(std);
+    }
+    else
+    {
+      bayesLinReg->Predict(testPoints, predictions);
+    }
+
+    // Save test predictions (one per line).
+    IO::GetParam<arma::mat>("predictions") = std::move(predictions);
+  }
+
+  IO::GetParam<BayesianLinearRegression*>("output_model") = bayesLinReg;
+}

src/mlpack/tests/CMakeLists.txt

@@ -9,6 +9,7 @@ add_executable(mlpack_test
   async_learning_test.cpp
   augmented_rnns_tasks_test.cpp
   binarize_test.cpp
+  bayesian_linear_regression_test.cpp
   callback_test.cpp
   cf_test.cpp
   cli_binding_test.cpp
@@ -100,6 +101,7 @@ add_executable(mlpack_test
   union_find_test.cpp
   vantage_point_tree_test.cpp
   wgan_test.cpp
+  main_tests/bayesian_linear_regression_test.cpp
   main_tests/cf_test.cpp
   main_tests/dbscan_test.cpp
   main_tests/decision_stump_test.cpp

src/mlpack/tests/bayesian_linear_regression_test.cpp

@@ -0,0 +1,194 @@
+/**
+ * @file tests/bayesian_linear_regression_test.cpp 
+ * @author Clement Mercier
+ *
+ * Test for BayesianLinearRegression.
+ *
+ * mlpack is free software; you may redistribute it and/or modify it under the
+ * terms of the 3-clause BSD license.  You should have received a copy of the
+ * 3-clause BSD license along with mlpack.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+
+#include <mlpack/core/data/load.hpp>
+#include <mlpack/methods/bayesian_linear_regression/bayesian_linear_regression.hpp>
+#include <mlpack/methods/linear_regression/linear_regression.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+using namespace mlpack::regression;
+using namespace mlpack::data;
+
+BOOST_AUTO_TEST_SUITE(BayesianLinearRegressionTest);
+
+void GenerateProblem(arma::mat& matX,
+                     arma::rowvec& y,
+                     size_t nPoints,
+                     size_t nDims,
+                     float sigma = 0.0)
+{
+  matX = arma::randn(nDims, nPoints);
+  arma::colvec omega = arma::randn(nDims);
+  // Compute y and add noise.
+  y = omega.t() * matX + arma::randn(nPoints).t() * sigma;
+}
+
+// Ensure that predictions are close enough to the target
+// for a free noise dataset.
+BOOST_AUTO_TEST_CASE(BayesianLinearRegressionRegressionTest)
+{
+  arma::mat matX;
+  arma::rowvec y, predictions;
+
+  GenerateProblem(matX, y, 200, 10);
+
+  // Instanciate and train the estimator.
+  BayesianLinearRegression estimator(true);
+  estimator.Train(matX, y);
+  estimator.Predict(matX, predictions);
+
+  // Check the predictions are close enough to the targets in a free noise case.
+  for (size_t i = 0; i < y.size(); i++)
+    BOOST_REQUIRE_CLOSE(predictions[i], y[i], 1e-6);
+
+  // Check that the estimated variance is zero.
+  BOOST_REQUIRE_SMALL(estimator.Variance(), 1e-6);
+}
+
+// Verify centerData and scaleData equal false do not affect the solution.
+BOOST_AUTO_TEST_CASE(TestCenter0ScaleData0)
+{
+  arma::mat matX;
+  arma::rowvec y;
+  size_t nDims = 30, nPoints = 100;
+
+  GenerateProblem(matX, y, nPoints, nDims, 0.5);
+
+  BayesianLinearRegression estimator(false, false);
+
+  estimator.Train(matX, y);
+
+  // Check dataOffset is empty.
+  BOOST_REQUIRE(estimator.DataOffset().n_elem == 0);
+
+  // To be neutral responseOffset must be 0.
+  BOOST_REQUIRE(estimator.ResponsesOffset() == 0);
+
+  // Check dataScale is empty.
+  BOOST_REQUIRE(estimator.DataScale().n_elem == 0);
+}
+
+// Verify that centering and normalization are correct.
+BOOST_AUTO_TEST_CASE(TestCenterDataTrueScaleDataTrue)
+{
+  arma::mat matX;
+  arma::rowvec y;
+  size_t nDims = 5, nPoints = 100;
+  GenerateProblem(matX, y, nPoints, nDims, 0.5);
+
+  BayesianLinearRegression estimator(true, true);
+  estimator.Train(matX, y);
+
+  arma::colvec xMean = arma::mean(matX, 1);
+  arma::colvec xStd = arma::stddev(matX, 0, 1);
+  double yMean = arma::mean(y);
+
+  BOOST_REQUIRE_SMALL((double) abs(sum(estimator.DataOffset() - xMean)), 1e-6);
+  BOOST_REQUIRE_SMALL((double) abs(sum(estimator.DataScale() - xStd)), 1e-6);
+  BOOST_REQUIRE_CLOSE(estimator.ResponsesOffset(), yMean, 1e-6);
+}
+
+// Make sure a model with center ans scale option set is different than a model
+// without it set.
+BOOST_AUTO_TEST_CASE(OptionsMakeModelDifferent)
+{
+  arma::mat matX;
+  arma::rowvec y;
+  size_t nDims = 10, nPoints = 100;
+  GenerateProblem(matX, y, nPoints, nDims, 0.5);
+
+  BayesianLinearRegression blr(false, false), blrC(true, false),
+      blrCS(true, true);
+
+  blr.Train(matX, y);
+  blrC.Train(matX, y);
+  blrCS.Train(matX, y);
+
+  for (size_t i = 0; i < nDims; ++i)
+    BOOST_REQUIRE((blr.Omega()(i) != blrC.Omega()(i)) &&
+                  (blr.Omega()(i) != blrCS.Omega()(i)) &&
+                  (blrC.Omega()(i) != blrCS.Omega()(i)));
+}
+
+// Check that Train() does not fail with two colinear vectors.
+BOOST_AUTO_TEST_CASE(SingularMatix)
+{
+  arma::mat matX;
+  arma::rowvec y;
+
+  GenerateProblem(matX, y, 200, 10);
+  // Now the first and the second rows are indentical.
+  matX.row(1) = matX.row(0);
+
+  BayesianLinearRegression estimator;
+  estimator.Train(matX, y);
+}
+
+// Check that std are well computed/coherent. At least higher than the
+// estimated predictive variance.
+BOOST_AUTO_TEST_CASE(PredictiveUncertainties)
+{
+  arma::mat matX;
+  arma::rowvec y;
+
+  GenerateProblem(matX, y, 100, 10, 1);
+
+  BayesianLinearRegression estimator(true, true);
+  estimator.Train(matX, y);
+
+  arma::rowvec responses, std;
+  estimator.Predict(matX, responses, std);
+  const double estStd = sqrt(estimator.Variance());
+
+  for (size_t i = 0; i < matX.n_cols; i++)
+    BOOST_REQUIRE_GT(std[i], estStd);
+
+  // Check that the estimated variance is close to 1.
+  BOOST_REQUIRE_CLOSE(estStd, 1, 30);
+}
+
+// Check the solution is equal to the classical ridge.
+BOOST_AUTO_TEST_CASE(EqualtoRidge)
+{
+  arma::mat matX;
+  arma::rowvec y, blrPred, ridgePred;
+
+  size_t trial = 0;
+  for ( ; trial < 3; ++trial)
+  {
+    GenerateProblem(matX, y, 100, 10, 1);
+
+    BayesianLinearRegression blr(false, false);
+    blr.Train(matX, y);
+
+    LinearRegression ridge(matX, y, blr.Alpha() / blr.Beta(), false);
+
+    blr.Predict(matX, blrPred);
+    ridge.Predict(matX, ridgePred);
+
+    // If the predictions seem far off, just try again.
+    if (arma::norm(blrPred - ridgePred) > 1e-5)
+      continue;
+
+    // Check the predictions are close enough between ridge and our blr.
+    for (size_t i = 0; i < y.size(); ++i)
+      BOOST_REQUIRE_CLOSE(blrPred[i], ridgePred[i], 1);
+
+    // Exit once a test case has completed.
+    break;
+  }
+
+  BOOST_REQUIRE_LT(trial, 3);
+}
+
+BOOST_AUTO_TEST_SUITE_END();

src/mlpack/tests/main_tests/bayesian_linear_regression_test.cpp

@@ -0,0 +1,143 @@
+/**
+ * @file tests/main_tests/bayesian_linear_regression_test.cpp
+ * @author Clement Mercier
+ *
+ * Test mlpackMain() of bayesian_linear_regression_main.cpp.
+ *
+ * mlpack is free software; you may redistribute it and/or modify it under the
+ * terms of the 3-clause BSD license.  You should have received a copy of the
+ * 3-clause BSD license along with mlpack.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#include <string>
+
+#define BINDING_TYPE BINDING_TYPE_TEST
+static const std::string testName = "BayesianLinearRegression";
+
+#include <mlpack/core.hpp>
+#include <mlpack/core/util/mlpack_main.hpp>
+#include "test_helper.hpp"
+#include <mlpack/methods/bayesian_linear_regression/bayesian_linear_regression_main.cpp>
+
+#include <boost/test/unit_test.hpp>
+#include "../test_tools.hpp"
+
+using namespace mlpack;
+
+struct BRTestFixture
+{
+ public:
+  BRTestFixture()
+  {
+    // Cache in the options for this program.
+    IO::RestoreSettings(testName);
+  }
+
+  ~BRTestFixture()
+  {
+    // Clear the settings.
+    bindings::tests::CleanMemory();
+    IO::ClearSettings();
+  }
+};
+
+BOOST_FIXTURE_TEST_SUITE(BayesianLinearRegressionMainTest, BRTestFixture);
+
+/**
+ * Check the center and scale options.
+ */
+BOOST_AUTO_TEST_CASE(BRCenter0Scale0)
+{
+  int n = 50, m = 4;
+  arma::mat matX = arma::randu<arma::mat>(m, n);
+  arma::rowvec omega = arma::randu<arma::rowvec>(m);
+  arma::rowvec y =  omega * matX;
+
+  SetInputParam("input", std::move(matX));
+  SetInputParam("responses", std::move(y));
+  SetInputParam("center", false);
+
+  mlpackMain();
+
+  BayesianLinearRegression* estimator =
+      IO::GetParam<BayesianLinearRegression*>("output_model");
+
+  BOOST_REQUIRE(estimator->DataOffset().n_elem == 0);
+  BOOST_REQUIRE(estimator->DataScale().n_elem == 0);
+}
+
+/**
+ * Check predictions of saved model and in code model are equal.
+ */
+BOOST_AUTO_TEST_CASE(BayesianLinearRegressionSavedEqualCode)
+{
+  int n = 10, m = 4;
+  arma::mat matX = arma::randu<arma::mat>(m, n);
+  arma::mat matXtest = arma::randu<arma::mat>(m, 2 * n);
+  const arma::rowvec omega = arma::randu<arma::rowvec>(m);
+  arma::rowvec y =  omega * matX;
+
+  BayesianLinearRegression model;
+  model.Train(matX, y);
+
+  arma::rowvec responses;
+  model.Predict(matXtest, responses);
+
+  SetInputParam("input", std::move(matX));
+  SetInputParam("responses", std::move(y));
+
+  mlpackMain();
+
+  IO::GetSingleton().Parameters()["input"].wasPassed = false;
+  IO::GetSingleton().Parameters()["responses"].wasPassed = false;
+
+  SetInputParam("input_model",
+                IO::GetParam<BayesianLinearRegression*>("output_model"));
+  SetInputParam("test", std::move(matXtest));
+
+  mlpackMain();
+
+  arma::mat ytest = std::move(responses);
+  // Check that initial output and output using saved model are same.
+  CheckMatrices(ytest, IO::GetParam<arma::mat>("predictions"));
+}
+
+/**
+ * Check a crash happens if neither input or input_model are specified.
+ * Check a crash happens if both input and input_model are specified.
+ */
+BOOST_AUTO_TEST_CASE(CheckParamsPassed)
+{
+  int n = 10, m = 4;
+  arma::mat matX = arma::randu<arma::mat>(m, n);
+  arma::mat matXtest = arma::randu<arma::mat>(m, 2 * n);
+  const arma::rowvec omega = arma::randu<arma::rowvec>(m);
+  arma::rowvec y =  omega * matX;
+
+  BayesianLinearRegression model;
+  model.Train(matX, y);
+
+  arma::rowvec responses;
+  model.Predict(matXtest, responses);
+
+  // Check that std::runtime_error is thrown if neither input or input_model
+  // is specified.
+  SetInputParam("responses", std::move(y));
+
+  BOOST_REQUIRE_THROW(mlpackMain(), std::runtime_error);
+
+  // Continue only with input passed.
+  SetInputParam("input", std::move(matX));
+  mlpackMain();
+
+  // Now pass the previous trained model and one input matrix at the same time.
+  // An error should occur.
+  SetInputParam("input", std::move(matX));
+  SetInputParam("input_model",
+                IO::GetParam<BayesianLinearRegression*>("output_model"));
+  SetInputParam("test", std::move(matXtest));
+
+  BOOST_REQUIRE_THROW(mlpackMain(), std::runtime_error);
+}
+
+BOOST_AUTO_TEST_SUITE_END();

src/mlpack/tests/serialization_test.cpp

@@ -39,6 +39,7 @@
 #include <mlpack/methods/lsh/lsh_search.hpp>
 #include <mlpack/methods/decision_stump/decision_stump.hpp>
 #include <mlpack/methods/lars/lars.hpp>
+#include <mlpack/methods/bayesian_linear_regression/bayesian_linear_regression.hpp>
 #include <mlpack/methods/ann/rbm/rbm.hpp>
 #include <mlpack/methods/ann/init_rules/gaussian_init.hpp>
 
@@ -1602,4 +1603,34 @@ BOOST_AUTO_TEST_CASE(ssRBMTest)
   CheckMatrices(Rbm.Weight(), RbmBinary.Weight());
 }
 
+// Make sure serialization works for BayesianLinearRegression.
+BOOST_AUTO_TEST_CASE(BayesianLinearRegressionTest)
+{
+  using namespace mlpack::regression;
+
+  // Create a dataset.
+  arma::mat matX = arma::randn(75, 250);
+  arma::vec omega = arma::randn(75, 1);
+  arma::rowvec y = omega.t() * matX;
+
+  BayesianLinearRegression blr(false, false);
+  blr.Train(matX, y);
+  arma::vec omegaOpt = blr.Omega();
+
+  // Now, serialize.
+  BayesianLinearRegression xmlBlr(false, false), binaryBlr(false, false),
+    textBlr(false, false);
+
+  SerializeObjectAll(blr, xmlBlr, binaryBlr, textBlr);
+
+  // Now, check that predictions are the same.
+  arma::rowvec pred, xmlPred, textPred, binaryPred;
+  blr.Predict(matX, pred);
+  xmlBlr.Predict(matX, xmlPred);
+  textBlr.Predict(matX, textPred);
+  binaryBlr.Predict(matX, binaryPred);
+
+  CheckMatrices(pred, xmlPred, textPred, binaryPred);
+}
+
 BOOST_AUTO_TEST_SUITE_END();