diff --git a/HISTORY.md b/HISTORY.md index a941b58491..67d39452d2 100644 --- a/HISTORY.md +++ b/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). diff --git a/src/mlpack/methods/CMakeLists.txt b/src/mlpack/methods/CMakeLists.txt index 83c96e68dd..d548d9c769 100644 --- a/src/mlpack/methods/CMakeLists.txt +++ b/src/mlpack/methods/CMakeLists.txt @@ -6,6 +6,7 @@ set(DIRS ann approx_kfn bias_svd + bayesian_linear_regression block_krylov_svd cf dbscan diff --git a/src/mlpack/methods/bayesian_linear_regression/CMakeLists.txt b/src/mlpack/methods/bayesian_linear_regression/CMakeLists.txt new file mode 100644 index 0000000000..5cdae4274d --- /dev/null +++ b/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") diff --git a/src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression.cpp b/src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression.cpp new file mode 100644 index 0000000000..10a5f92bd5 --- /dev/null +++ b/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 +#include + +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(data.memptr()), data.n_rows, + data.n_cols, false, true); + responsesProc = arma::rowvec(const_cast(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(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(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; + } +} diff --git a/src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression.hpp b/src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression.hpp new file mode 100644 index 0000000000..415f63e595 --- /dev/null +++ b/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 + +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 + 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 diff --git a/src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression_impl.hpp b/src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression_impl.hpp new file mode 100644 index 0000000000..9c4cb20f09 --- /dev/null +++ b/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 +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 diff --git a/src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression_main.cpp b/src/mlpack/methods/bayesian_linear_regression/bayesian_linear_regression_main.cpp new file mode 100644 index 0000000000..86fc946fe2 --- /dev/null +++ b/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 +#include +#include + +#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("center"); + bool scale = IO::GetParam("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("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("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("input_model"); + } + + if (IO::HasParam("test")) + { + Log::Info << "Regressing on test points." << endl; + // Load test points. + mat testPoints = std::move(IO::GetParam("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("stds") = std::move(std); + } + else + { + bayesLinReg->Predict(testPoints, predictions); + } + + // Save test predictions (one per line). + IO::GetParam("predictions") = std::move(predictions); + } + + IO::GetParam("output_model") = bayesLinReg; +} diff --git a/src/mlpack/tests/CMakeLists.txt b/src/mlpack/tests/CMakeLists.txt index 277342fb00..e1a6768c04 100644 --- a/src/mlpack/tests/CMakeLists.txt +++ b/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 diff --git a/src/mlpack/tests/bayesian_linear_regression_test.cpp b/src/mlpack/tests/bayesian_linear_regression_test.cpp new file mode 100644 index 0000000000..d0c34509b9 --- /dev/null +++ b/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 +#include +#include + +#include + +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(); diff --git a/src/mlpack/tests/main_tests/bayesian_linear_regression_test.cpp b/src/mlpack/tests/main_tests/bayesian_linear_regression_test.cpp new file mode 100644 index 0000000000..e1b6773a5e --- /dev/null +++ b/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 + +#define BINDING_TYPE BINDING_TYPE_TEST +static const std::string testName = "BayesianLinearRegression"; + +#include +#include +#include "test_helper.hpp" +#include + +#include +#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(m, n); + arma::rowvec omega = arma::randu(m); + arma::rowvec y = omega * matX; + + SetInputParam("input", std::move(matX)); + SetInputParam("responses", std::move(y)); + SetInputParam("center", false); + + mlpackMain(); + + BayesianLinearRegression* estimator = + IO::GetParam("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(m, n); + arma::mat matXtest = arma::randu(m, 2 * n); + const arma::rowvec omega = arma::randu(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("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("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(m, n); + arma::mat matXtest = arma::randu(m, 2 * n); + const arma::rowvec omega = arma::randu(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("output_model")); + SetInputParam("test", std::move(matXtest)); + + BOOST_REQUIRE_THROW(mlpackMain(), std::runtime_error); +} + +BOOST_AUTO_TEST_SUITE_END(); diff --git a/src/mlpack/tests/serialization_test.cpp b/src/mlpack/tests/serialization_test.cpp index e30fcc4dc9..327786e81c 100644 --- a/src/mlpack/tests/serialization_test.cpp +++ b/src/mlpack/tests/serialization_test.cpp @@ -39,6 +39,7 @@ #include #include #include +#include #include #include @@ -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();