apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'upstream/master' into batch

Browse Source
This commit is contained in:
Ryan Curtin
2017-10-09 14:52:41 -04:00
parent e5243262fe 8b5765bb35
commit e6ecfb897c
19 changed files with 1492 additions and 271 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/guide/cv.hpp
+372
View File
@@ -0,0 +1,372 @@
namespace mlpack {
namespace cv {
/*! @page cv Cross-Validation
@section cvintro Introduction
@b mlpack implements cross-validation support for its learning algorithms, for a
variety of performance measures. Cross-validation is useful for determining an
estimate of how well the learner will generalize to un-seen test data. It is a
commonly used part of the data science pipeline.
In short, given some learner and some performance measure, we wish to get an
average of the performance measure given different splits of the dataset into
training data and validation data. The learner is trained on the training data,
and the performance measure is evaluated on the validation data.
mlpack currently implements two easy-to-use forms of cross-validation:
- @b simple @b cross-validation, where we simply desire the performance measure
on a single split of the data into a training set and validation set
- @b k-fold @b cross-validation, where we split the data k ways and desire the
average performance measure on each of the k splits of the data
In this tutorial we will see the usage examples and details of the
cross-validation module. Because the cross-validation code is generic and can
be used with any learner and performance measure, any use of the
cross-validation code in mlpack has to be in C++.
This tutorial is split into the following sections:
- @ref cvbasic Simple cross-validation examples
- @ref cvbasic_ex_1 10-fold cross-validation on softmax regression
- @ref cvbasic_ex_2 10-fold cross-validation on weighted decision trees
- @ref cvbasic_ex_3 10-fold cross-validation with categorical decision trees
- @ref cvbasic_ex_4 Simple cross-validation for linear regression
- @ref cvbasic_metrics Performance measures
- @ref cvbasic_api The \c KFoldCV and \c SimpleCV classes
- @ref cvbasic_further Further reference
@section cvbasic Simple cross-validation examples
@subsection cvbasic_ex_1 10-fold cross-validation on softmax regression
Suppose we have some data to train and validate on, as defined below:
@code
// 100-point 6-dimensional random dataset.
arma::mat data = arma::randu<arma::mat>(6, 100);
// Random labels in the [0, 4] interval.
arma::Row<size_t> labels =
arma::randi<arma::Row<size_t>>(100, arma::distr_param(0, 4));
size_t numClasses = 5;
@endcode
The code above generates an 100-point random 6-dimensional dataset with 5
classes.
To run 10-fold cross-validation for softmax regression with accuracy as a
performance measure, we can write the following piece of code.
@code
KFoldCV<SoftmaxRegression, Accuracy> cv(10, data, labels, numClasses);
double lambda = 0.1;
double softmaxAccuracy = cv.Evaluate(lambda);
@endcode
Note that the \c Evaluate method of \c KFoldCV takes any hyperparameters of an
algorithm---that is, anything that is not \c data, \c labels, \c numClasses,
\c datasetInfo, or \c weights (those last three may not be present for every
algorithm type). To be more specific, in this example the \c Evaluate method
relies on the following \ref regression::SoftmaxRegression "SoftmaxRegression"
constructor:
@code
template<typename OptimizerType = mlpack::optimization::L_BFGS>
SoftmaxRegression(const arma::mat& data,
const arma::Row<size_t>& labels,
const size_t numClasses,
const double lambda = 0.0001,
const bool fitIntercept = false,
OptimizerType optimizer = OptimizerType());
@endcode
which has the parameter \c lambda after three conventional arguments (\c data,
\c labels and \c numClasses). We can skip passing \c fitIntercept and \c
optimizer since there are the default values. (Technically, we don't even need
to pass \c lambda since there is a default value.)
In general to cross-validate you need to specify what machine learning algorithm
and metric you are going to use, and then to pass some conventional data-related
parameters into one of the cross-validation constructors and all other
parameters (which are generally hyperparameters) into the \c Evaluate method.
@subsection cvbasic_ex_2 10-fold cross-validation on weighted decision trees
In the following example we will cross-validate
\ref tree::DecisionTree "DecisionTree" with weights. This is very similar to
the previous example, except that we also have instance weights for each point
in the dataset. We can generate weights for the dataset from the previous
example with the code below:
@code
// Random weights for every point from the code snippet above.
arma::rowvec weights = arma::randu<arma::mat>(1, 100);
@endcode
Given those weights for each point, we can now perform cross-validation by also
passing the weights to the constructor of \c KFoldCV:
@code
KFoldCV<DecisionTree<>, Accuracy> cv2(10, data, labels, numClasses, weights);
size_t minimumLeafSize = 8;
double weightedDecisionTreeAccuracy = cv2.Evaluate(minimumLeafSize);
@endcode
As with the previous example, internally this call to \c cv2.Evaluate() relies
on the following \ref tree::DecisionTree "DecisionTree" constructor:
@code
template<typename MatType, typename LabelsType, typename WeightsType>
DecisionTree(MatType&& data,
LabelsType&& labels,
const size_t numClasses,
WeightsType&& weights,
const size_t minimumLeafSize = 10,
const std::enable_if_t<arma::is_arma_type<
typename std::remove_reference<WeightsType>::type>::value>*
= 0);
@endcode
@subsection cvbasic_ex_3 10-fold cross-validation with categorical decision trees
\ref tree::DecisionTree "DecisionTree" models can be constructed in multiple
other ways. For example, if we have a dataset with both categorical and
numerical features, we can also perform cross-validation by using the associated
\c data::DatasetInfo object. Thus, given some \c data::DatasetInfo object
called \c datasetInfo (that perhaps was produced by a call to \c data::Load() ),
we can perform k-fold cross-validation in a similar manner to the other
examples:
@code
KFoldCV<DecisionTree<>, Accuracy> cv3(10, data, datasetInfo, labels,
numClasses);
double decisionTreeWithDIAccuracy = cv3.Evaluate(minimumLeafSize);
@endcode
This particular call to \c cv3.Evaluate() relies on the following
\ref tree::DecisionTree "DecisionTree" constructor:
@code
template<typename MatType, typename LabelsType>
DecisionTree(MatType&& data,
const data::DatasetInfo& datasetInfo,
LabelsType&& labels,
const size_t numClasses,
const size_t minimumLeafSize = 10);
@endcode
@subsection cvbasic_ex_4 Simple cross-validation for linear regression
\c SimpleCV has the same interface as \c KFoldCV, except it takes as one of its
arguments a proportion (from 0 to 1) of data used as a validation set. For
example, to validate \ref regression::LinearRegression "LinearRegression" with
20\% of the data used in the validation set we can write the following code.
@code
// Random responses for every point from the code snippet in the beginning of
// the tutorial.
arma::rowvec responses = arma::randu<arma::rowvec>(100);
SimpleCV<LinearRegression, MSE> cv4(0.2, data, responses);
double lrLambda = 0.05;
double lrMSE = cv4.Evaluate(lrLambda);
@endcode
@section cvbasic_metrics Performance measures
The cross-validation classes require a performance measure to be specified.
\b mlpack has a number of performance measures implemented; below is a list:
- mlpack::cv::Accuracy: a simple measure of accuracy
- mlpack::cv::F1: the F1 score; depends on an averaging strategy
- mlpack::cv::MSE: minimum squared error (for regression problems)
- mlpack::cv::Precision: the precision, for classification problems
- mlpack::cv::Recall: the recall, for classification problems
In addition, it is not difficult to implement a custom performance measure. A
class following the structure below can be used:
@code
class CustomMeasure
{
//
// This evaluates the metric given a trained model and a set of data (with
// labels or responses) to evaluate on. The data parameter will be a type of
// Armadillo matrix, and the labels will be the labels that go with the model.
//
// If you know that your model is a classification model (and thus that
// ResponsesType will be arma::Row<size_t>), it is ok to replace the
// ResponsesType template parameter with arma::Row<size_t>.
//
template<typename MLAlgorithm, typename DataType, typename ResponsesType>
static double Evaluate(MLAlgorithm& model,
const DataType& data,
const ResponsesType& labels)
{
// Inside the method you should call model.Predict() and compare the
// values with the labels, in order to get the desired performance measure
// and return it.
}
};
@endcode
Once this is implemented, then \c CustomMeasure (or whatever the class is
called) is easy to use as a custom performance measure with \c KFoldCV or
\c SimpleCV.
@section cvbasic_api The KFoldCV and SimpleCV classes
This section provides details about the \c KFoldCV and \c SimpleCV classes.
The cross-validation infrastructure is based on heavy amounts of template
metaprogramming, so that any \b mlpack learner and any performance measure can
be used. Both classes have two required template parameters and one optional
parameter:
- \c MLAlgorithm: the type of learner to be used
- \c Metric: the performance measure to be evaluated
- \c MatType: the type of matrix used to store the data
In addition, there are two more template parameters, but these are automatically
extracted from the given \c MLAlgorithm class, and users should not need to
specify these parameters except when using an unconventional type like
\c arma::fmat for data points.
The general structure of the \c KFoldCV and \c SimpleCV classes is split into
two parts:
- The constructor: create the object, and store the data for the \c MLAlgorithm
training.
- The \c Evaluate() method: take any non-data parameters for the
\c MLAlgorithm and calculate the desired performance measure.
This split is important because it defines the API: all data-related parameters
are passed to the constructor, whereas algorithm hyperparameters are passed to
the \c Evaluate() method.
@subsection cvbasic_api_constructor The KFoldCV and SimpleCV constructors
There are six constructors available for \c KFoldCV and \c SimpleCV, each
tailored for a different learning situation. Each is given below for the
\c KFoldCV class, but the same constructors are also available for the
\c SimpleCV class, with the exception that instead of specifying \c k, the
number of folds, the \c SimpleCV class takes a parameter between 0 and 1
specifying the percentage of the dataset to use as a validation set.
- `KFoldCV(k, xs, ys)`: this is for unweighted regression applications and
two-class classification applications; \c xs is the dataset and \c ys
are the responses or labels for each point in the dataset.
- `KFoldCV(k, xs, ys, numClasses)`: this is for unweighted classification
applications; \c xs is the dataset, \c ys are the class labels for each
data point, and \c numClasses is the number of classes in the dataset.
- `KFoldCV(k, xs, datasetInfo, ys, numClasses)`: this is for unweighted
categorical/numeric classification applications; \c xs is the dataset,
\c datasetInfo is a data::DatasetInfo object that holds the types of
each dimension in the dataset, \c ys are the class labels for each data
point, and \c numClasses is the number of classes in the dataset.
- `KFoldCV(k, xs, ys, weights)`: this is for weighted regression or
two-class classification applications; \c xs is the dataset, \c ys are
the responses or labels for each point in the dataset, and \c weights
are the weights for each point in the dataset.
- `KFoldCV(k, xs, ys, numClasses, weights)`: this is for weighted
classification applications; \c xs is the dataset, \c ys are the class
labels for each point in the dataset; \c numClasses is the number of
classes in the dataset, and \c weights holds the weights for each point
in the dataset.
- `KFoldCV(k, xs, datasetInfo, ys, numClasses, weights)`: this is for
weighted cateogrical/numeric classification applications; \c xs is the
dataset, \c datasetInfo is a data::DatasetInfo object that holds the
types of each dimension in the dataset, \c ys are the class labels for
each data point, \c numClasses is the number of classes in each dataset,
and \c weights holds the weights for each point in the dataset.
Note that the constructor you should use is the constructor that most closely
matches the constructor of the machine learning algorithm you would like
performance measures of. So, for instance, if you are doing multi-class softmax
regression, you could call the constructor
\c "SoftmaxRegression(xs, ys, numClasses)". Therefore, for \c KFoldCV you would
call the constructor \c "KFoldCV(k, xs, ys, numClasses)" and for \c SimpleCV you
would call the constructor \c "SimpleCV(pct, xs, ys, numClasses)".
@subsection cvbasic_api_evaluate The Evaluate() method
The other method that \c KFoldCV and \c SimpleCV have is the method to
actually calculate the performance measure: \c Evaluate(). The \c Evaluate()
method takes any hyperparameters that would follow the data arguments to the
constructor or \c Train() method of the given \c MLAlgorithm. The
\c Evaluate() method takes no more arguments than that, and returns the
desired performance measure on the dataset.
Therefore, let us suppose that we are interested in cross-validating the
performance of a softmax regression model, and that we have constructed
the appropriate \c KFoldCV object using the code below:
@code
KFoldCV<SoftmaxRegression, Precision> cv(k, data, labels, numClasses);
@endcode
The \ref regression::SoftmaxRegression "SoftmaxRegression" class has the
constructor
@code
template<typename OptimizerType = mlpack::optimization::L_BFGS>
SoftmaxRegression(const arma::mat& data,
const arma::Row<size_t>& labels,
const size_t numClasses,
const double lambda = 0.0001,
const bool fitIntercept = false,
OptimizerType optimizer = OptimizerType());
@endcode
Note that all parameters after are \c numClasses are optional. This means that
we can specify none or any of them in our call to \c Evaluate(). Below is some
example code showing three different ways we can call \c Evaluate() with the
\c cv object from the code snippet above.
@code
// First, call with all defaults.
double result1 = cv.Evaluate();
// Next, call with lambda set to 0.1 and fitIntercept set to true.
double result2 = cv.Evaluate(0.1, true);
// Lastly, create a custom optimizer to use for optimization, and use a lambda
// value of 0.5 and fit no intercept.
optimization::SGD<> sgd(0.05, 50000); // Step size of 0.05, 50k max iterations.
double result3 = cv.Evaluate(0.5, false, sgd);
@endcode
The same general idea applies to any \c MLAlgorithm: all hyperparameters must be
passed to the \c Evaluate() method of \c KFoldCV or \c SimpleCV.
@section cvbasic_further Further references
For further documentation, please see the associated Doxygen documentation for
each of the relevant classes:
- mlpack::cv::SimpleCV
- mlpack::cv::KFoldCV
- mlpack::cv::Accuracy
- mlpack::cv::F1
- mlpack::cv::MSE
- mlpack::cv::Precision
- mlpack::cv::Recall
If you are interested in implementing a different cross-validation strategy than
k-fold cross-validation or simple cross-validation, take a look at the
implementations of each of those classes to guide your implementation.
In addition, the @ref hpt "hyperparameter tuner" documentation may also be
relevant.
*/
} // namespace cv
} // namespace mlpack
doc/guide/hpt.hpp
+238
View File
@@ -0,0 +1,238 @@
namespace mlpack {
namespace hpt {
/*! @page hpt Hyper-Parameter Tuning
@section hptintro Introduction
\b mlpack implements a generic hyperparameter tuner that is able to tune both
continuous and discrete parameters of various different algorithms. This is an
important task---the performance of many machine learning algorithms can be
highly dependent on the hyperparameters that are chosen for that algorithm.
(One example: the choice of \f$k\f$ for a \f$k\f$-nearest-neighbors classifier.)
This hyper-parameter tuner is built on the same general concept as the
cross-validation classes (see the @ref cv "cross-validation tutorial"): given
some machine learning algorithm, some data, some performance measure, and a set
of hyperparameters, attempt to find the hyperparameter set that best optimizes
the performance measure on the given data with the given algorithm.
\b mlpack's implementation of hyperparameter tuning is flexible, and is built in
a way that supports many algorithms and many optimizers. At the time of this
writing, complex hyperparameter optimization techniques are not available, but
the hyperparameter tuner does support these, should they be implemented in the
future.
In this tutorial we will see the usage examples of the hyper-parameter tuning
module, and also more details about the \c HyperParameterTuner class.
@section hptbasic Basic Usage
The interface of the hyper-parameter tuning module is quite similar to the
interface of the @ref cv "cross-validation module". To construct a \c
HyperParameterTuner object you need to specify as template parameters what
machine learning algorithm, cross-validation strategy, performance measure, and
optimization strategy (\ref optimization::GridSearch "GridSearch" will be used by
default) you are going to use. Then, you must pass the same arguments as for
the cross-validation classes: the data and labels (or responses) to use are
given to the constructor, and the possible hyperparameter values are given to
the \c HyperParameterTuner::Optimize() method, which returns the best
algorithm configuration as a \c std::tuple<>.
Let's see some examples.
Suppose we have the following data to train and validate on.
@code
// 100-point 5-dimensional random dataset.
arma::mat data = arma::randu<arma::mat>(5, 100);
// Noisy responses retrieved by a random linear transformation of data.
arma::rowvec responses = arma::randu<arma::rowvec>(5) * data +
0.1 * arma::randn<arma::rowvec>(100);
@endcode
Given the dataset above, we can use the following code to try to find a good \c
lambda value for \ref regression::LinearRegression "LinearRegression". Here we
use \ref cv::SimpleCV "SimpleCV" instead of k-fold cross-validation to save
computation time.
@code
// Using 80% of data for training and remaining 20% for assessing MSE.
double validationSize = 0.2;
HyperParameterTuner<LinearRegression, MSE, SimpleCV> hpt(validationSize,
data, responses);
// Finding a good value for lambda from the discrete set of values 0.0, 0.001,
// 0.01, 0.1, and 1.0.
arma::vec lambdas{0.0, 0.001, 0.01, 0.1, 1.0};
double bestLambda;
std::tie(bestLambda) = hpt.Optimize(lambdas);
@endcode
In this example we have used \ref optimization::GridSearch "GridSearch" (the
default optimizer) to find a good value for the \c lambda hyper-parameter. For
that we have specified what values should be tried.
@section hptfixed Fixed Arguments
When some hyper-parameters should not be optimized, you can specify values
for them with the \c Fixed() method as in the following example of trying to
find good \c lambda1 and \c lambda2 values for \ref regression::LARS "LARS"
(least-angle regression).
@code
HyperParameterTuner<LARS, MSE, SimpleCV> hpt2(validationSize, data,
responses);
// The hyper-parameter tuner should not try to change the transposeData or
// useCholesky parameters.
bool transposeData = true;
bool useCholesky = false;
// We wish only to search for the best lambda1 and lambda2 values.
arma::vec lambda1Set{0.0, 0.001, 0.01, 0.1, 1.0};
arma::vec lambda2Set{0.0, 0.002, 0.02, 0.2, 2.0};
double bestLambda1, bestLambda2;
std::tie(bestLambda1, bestLambda2) = hpt2.Optimize(Fixed(transposeData),
Fixed(useCholesky), lambda1Set, lambda2Set);
@endcode
Note that for the call to \c hpt2.Optimize(), we have used the same order of
arguments as they appear in the corresponding \ref regression::LARS "LARS"
constructor:
@code
LARS(const arma::mat& data,
const arma::rowvec& responses,
const bool transposeData = true,
const bool useCholesky = false,
const double lambda1 = 0.0,
const double lambda2 = 0.0,
const double tolerance = 1e-16);
@endcode
@section hptgradient Gradient-Based Optimization
In some cases we may wish to optimize a hyperparameter over the space of all
possible real values, instead of providing a grid in which to search.
Alternately, we may know approximately optimal values from a grid search for
real-valued hyperparameters, but wish to further tune those values.
In this case, we can use a gradient-based optimizer for hyperparameter search.
In the following example, we try to optimize the \c lambda1 and \c lambda2
hyper-parameters for \ref regression::LARS "LARS" with the
\ref optimization::GradientDescent "GradientDescent" optimizer.
@code
HyperParameterTuner<LARS, MSE, SimpleCV, GradientDescent> hpt3(validationSize,
data, responses);
// GradientDescent can be adjusted in the following way.
hpt3.Optimizer().StepSize() = 0.1;
hpt3.Optimizer().Tolerance() = 1e-15;
// We can set up values used for calculating gradients.
hpt3.RelativeDelta() = 0.01;
hpt3.MinDelta() = 1e-10;
double initialLambda1 = 0.001;
double initialLambda2 = 0.002;
double bestGDLambda1, bestGDLambda2;
std::tie(bestGDLambda1, bestGDLambda2) = hpt3.Optimize(Fixed(transposeData),
Fixed(useCholesky), initialLambda1, initialLambda2);
@endcode
@section hpt_class The HyperParameterTuner class
The \c HyperParameterTuner class is very similar to the
\ref cv::KFoldCV "KFoldCV" and \ref cv::SimpleCV "SimpleCV" classes (see the
@ref "cross-validation tutorial" for more information on those two classes), but
there are a few important differences.
First, the \c HyperParameterTuner accepts five different hyperparameters; only
the first three of these are required:
- \c MLAlgorithm This is the algorithm to be used.
- \c Metric This is the performance measure to be used; see
@ref cvbasic_metrics for more information.
- \c CVType This is the type of cross-validation to be used for evaluating the
performance measure; this should be \ref cv::KFoldCV "KFoldCV" or
\ref cv::SimpleCV "SimpleCV".
- \c OptimizerType This is the type of optimizer to use; it can be
\c GridSearch or a gradient-based optimizer.
- \c MatType This is the type of data matrix to use. The default is
\c arma::mat. This only needs to be changed if you are specifically
using sparse data, or if you want to use a numeric type other than
\c double.
The last two template parameters are automatically inferred by the
\c HyperParameterTuner and should not need to be manually specified, unless an
unconventional data type like \c arma::fmat is being used for data points.
Typically, \ref cv::SimpleCV "SimpleCV" is a good choice for \c CVType because
it takes so much less time to compute than full \ref cv::KFoldCV "KFoldCV";
however, the disadvantage is that \ref cv::SimpleCV "SimpleCV" might give a
somewhat more noisy estimate of the performance measure on unseen test data.
The constructor for the \c HyperParameterTuner is called with exactly the same
arguments as the corresponding \c CVType that has been chosen. For more
information on that, please see the
@ref cvbasic_api "cross-validation constructor tutorial". As an example, if we
are using \ref cv::SimpleCV "SimpleCV" and wish to hold out 20\% of the dataset
as a validation set, we might construct a \c HyperParameterTuner like this:
@code
// We will use LinearRegression as the MLAlgorithm, and MSE as the performance
// measure. Our dataset is 'dataset' and the responses are 'responses'.
HyperParameterTuner<LinearRegression, MSE, SimpleCV> hpt(0.2, dataset,
responses);
@endcode
Next, we must set up the hyperparameters to be optimized. If we are doing a
grid search with the \ref optimization::GridSearch "GridSearch" optimizer (the
default), then we only need to pass a `std::vector` (for non-numeric
hyperparameters) or an `arma::vec` (for numeric hyperparameters) containing all
of the possible choices that we wish to search over.
For instance, a set of numeric values might be chosen like this, for the
\c lambda parameter (of type \c double):
@code
arma::vec lambdaSet = arma::vec("0.0 0.1 0.5 1.0");
@endcode
Similarly, a set of non-numeric values might be chosen like this, for the
\c intercept parameter:
@code
std::vector<bool> interceptSet = { false, true };
@endcode
Once all of these are set up, the \c HyperParameterTuner::Optimize() method may
be called to find the best set of hyperparameters:
@code
bool intercept;
double lambda;
std::tie(lambda, intercept) = hpt.Optimize(lambdaSet, interceptSet);
@endcode
Alternately, the \c Fixed() method (detailed in the @ref hptfixed
"Fixed arguments" section) can be used to fix the values of some parameters.
For continuous optimizers like
\ref optimization::GradientDescent "GradientDescent", a range does not need to
be specified but instead only a single value. See the
\ref hptgradient "Gradient-Based Optimization" section for more details.
@section hptfurther Further documentation
For more information on the \c HyperParameterTuner class, see the
mlpack::hpt::HyperParameterTuner class documentation and the
@ref cv "cross-validation tutorial".
*/
} // namespace hpt
} // namespace mlpack
doc/tutorials/tutorials.txt
+2
View File
@@ -18,6 +18,8 @@ start.
- \ref iodoc
- \ref timer
- \ref sample
- \ref cv
- \ref hpt
@section method_tut Method-specific Tutorials
src/mlpack/core.hpp
+2
View File
@@ -100,6 +100,8 @@
* - @ref iodoc
* - @ref timer
* - @ref sample
* - @ref cv
* - @ref hpt
* - @ref verinfo
*
* Tutorials on specific methods are also available.
src/mlpack/core/tree/CMakeLists.txt
+1
View File
@@ -120,6 +120,7 @@ set(SOURCES
statistic.hpp
traversal_info.hpp
tree_traits.hpp
enumerate_tree.hpp
)
# add directory name to sources
src/mlpack/core/tree/cosine_tree/cosine_tree.cpp
+2 -1
View File
@@ -94,7 +94,8 @@ CosineTree::CosineTree(const arma::mat& dataset,
// Initialize Monte Carlo error estimate for comparison.
double monteCarloError = root.FrobNormSquared();
while (monteCarloError > epsilon * root.FrobNormSquared())
while (treeQueue.top() &&
(monteCarloError > epsilon * root.FrobNormSquared()))
{
// Pop node from queue with highest projection error.
CosineTree* currentNode;
src/mlpack/core/tree/enumerate_tree.hpp
+64
View File
@@ -0,0 +1,64 @@
/**
* @file enumerate_tree.hpp
* @author Ivan (Jonan) Georgiev
*
* This file contains function that performs a simple depth-first walk on the tree
* calling `Enter` and `Leave` methods of a provided walker.
*
* 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_CORE_TREE_ENUMERATE_TREE_HPP
#define MLPACK_CORE_TREE_ENUMERATE_TREE_HPP
namespace mlpack {
namespace tree /** Trees and tree-building procedures. */ {
namespace enumerate {
// Actual implementation of the enumeration. The problem is the unified
// detection if we're on the root, because Enter and Leave expect the
// parent being passed.
template <class TreeType, class Walker>
void EnumerateTreeImpl(TreeType* tree, Walker& walker, bool root)
{
if (root)
walker.Enter(tree, (const TreeType*)nullptr);
const size_t numChildren = tree->NumChildren();
for (size_t i = 0; i < numChildren; ++i)
{
TreeType* child = tree->ChildPtr(i);
walker.Enter(child, tree);
EnumerateTreeImpl(child, walker, false);
walker.Leave(child, tree);
}
if (root)
walker.Leave(tree, (const TreeType*)nullptr);
}
} // namespace enumerate
/**
* Traverses all nodes of the tree, including the inner ones. On each node
* two methods of the `enumer` are called:
*
* Enter(TreeType* node, TreeType* parent);
* Leave(TreeType* node, TreeType* parent);
*
* @param walker An instance of custom class, receiver of the enumeration.
*/
template <class TreeType, class Walker>
inline void EnumerateTree(TreeType* tree, Walker& walker)
{
enumerate::EnumerateTreeImpl(tree, walker, true);
}
} // namespace tree
} // namespace mlpack
#endif // MLPACK_CORE_TREE_ENUMERATE_TREE_HPP
src/mlpack/core/util/timers.cpp
-1
View File
@@ -111,7 +111,6 @@ void Timers::PrintTimer(const std::string& timerName)
Log::Info << ", ";
Log::Info << s.count() << "." << std::setw(1)
<< (totalDurationMicroSec.count() / 100000) << " secs";
output = true;
}
Log::Info << ")";
src/mlpack/methods/det/det_main.cpp
+162 -44
View File
@@ -33,6 +33,15 @@ PROGRAM_INFO("Density Estimation With Density Estimation Trees",
" with the " + PRINT_PARAM_STRING("training_set_estimates") + " output "
"parameter."
"\n\n"
"Enabling path printing for each node outputs the path from the root node "
"to a leaf for each entry in the test set, or training set (if a test set "
"is not provided). Strings like 'LRLRLR' (indicating that traversal went "
"to the left child, then the right child, then the left child, and so "
"forth) will be output. If 'lr-id' or 'id-lr' are given as the " +
PRINT_PARAM_STRING("path_format") + " parameter, then the ID (tag) of "
"every node along the path will be printed after or before the L or R "
"character indicating the direction of traversal, respectively."
"\n\n"
"This program also can provide density estimates for a set of test points, "
"specified in the " + PRINT_PARAM_STRING("test") + " parameter. The "
"density estimation tree used for this task will be the tree that was "
@@ -61,6 +70,19 @@ PARAM_MATRIX_OUT("test_set_estimates", "The output estimates on the test set "
PARAM_MATRIX_OUT("vi", "The output variable importance values for each "
"feature.", "i");
// Tagging and path printing options
PARAM_STRING_IN("path_format", "The format of path printing: 'lr', 'id-lr', or "
"'lr-id'.", "p", "lr");
PARAM_STRING_OUT("tag_counters_file", "The file to output the number of points "
"that went to each leaf.", "c");
PARAM_STRING_OUT("tag_file", "The file to output the tags (and possibly paths)"
" for each sample in the test set.", "g");
PARAM_FLAG("skip_pruning", "Whether to bypass the pruning process and output "
"the unpruned tree only.", "s");
// Parameters for the training algorithm.
PARAM_INT_IN("folds", "The number of folds of cross-validation to perform for "
"the estimation (0 is LOOCV)", "f", 10);
@@ -76,71 +98,80 @@ PARAM_FLAG("volume_regularization", "This flag gives the used the option to use"
"penalize low volume leaves.", "R");
*/
void mlpackMain()
{
// Validate input parameters.
if (CLI::HasParam("training") && CLI::HasParam("input_model"))
Log::Fatal << "Only one of --training_file (-t) or --input_model_file (-m) "
<< "may be specified!" << endl;
Log::Fatal << "Only one of " << PRINT_PARAM_STRING("training") << " or " <<
PRINT_PARAM_STRING("input_model") << " may be specified!" << endl;
if (!CLI::HasParam("training") && !CLI::HasParam("input_model"))
Log::Fatal << "Neither --training_file (-t) nor --input_model_file (-m) "
<< "are specified!" << endl;
Log::Fatal << "Neither " << PRINT_PARAM_STRING("training") << " nor " <<
PRINT_PARAM_STRING("input_model") << " are specified!" << endl;
if (CLI::HasParam("tag_file") &&
!CLI::HasParam("training") && !CLI::HasParam("test"))
{
Log::Fatal << "Neither " << PRINT_PARAM_STRING("training") << " nor " <<
PRINT_PARAM_STRING("test") << " are specified, but needed when " <<
PRINT_PARAM_STRING("tag_file") << " is asked." << endl;
}
if (!CLI::HasParam("training"))
{
if (CLI::HasParam("training_set_estimates"))
Log::Warn << "--training_set_estimates_file (-e) ignored because "
<< "--training_file (-t) is not specified." << endl;
Log::Warn << PRINT_PARAM_STRING("training_set_estimates") <<
" ignored because " << PRINT_PARAM_STRING("training") <<
" is not specified." << endl;
if (CLI::HasParam("folds"))
Log::Warn << "--folds (-f) ignored because --training_file (-t) is not "
<< "specified." << endl;
Log::Warn << PRINT_PARAM_STRING("folds") << " ignored because " <<
PRINT_PARAM_STRING("training") << " is not specified." << endl;
if (CLI::HasParam("min_leaf_size"))
Log::Warn << "--min_leaf_size (-l) ignored because --training_file (-t) "
<< "is not specified." << endl;
Log::Warn << PRINT_PARAM_STRING("min_leaf_size") << " ignored because " <<
PRINT_PARAM_STRING("training") << " is not specified." << endl;
if (CLI::HasParam("max_leaf_size"))
Log::Warn << "--max_leaf_size (-L) ignored because --training_file (-t) "
<< "is not specified." << endl;
Log::Warn << PRINT_PARAM_STRING("max_leaf_size") << " ignored because " <<
PRINT_PARAM_STRING("training") << " is not specified." << endl;
}
else if (!CLI::HasParam("output_model") &&
!CLI::HasParam("training_set_estimates") &&
!CLI::HasParam("vi"))
{
Log::Warn << "None of --output_model_file (-M), --training_set_estimates "
<< "(-e), or --vi (-i) are specified; no output will be saved!" << endl;
Log::Warn << "None of " << PRINT_PARAM_STRING("output_model") << ", " <<
PRINT_PARAM_STRING("training_set_estimates") << ", or " <<
PRINT_PARAM_STRING("vi") << " are specified; no output will be saved!" <<
endl;
}
if (!CLI::HasParam("test") && CLI::HasParam("test_set_estimates"))
Log::Warn << "--test_set_estimates_file (-E) ignored because --test_file "
<< "(-T) is not specified." << endl;
Log::Warn << PRINT_PARAM_STRING("test_set_estimates") << " ignored " <<
"because " << PRINT_PARAM_STRING("test") << " is not specified." << endl;
// Are we training a DET or loading from file?
DTree<arma::mat, int>* tree;
arma::mat trainingData;
arma::mat testData;
if (CLI::HasParam("training"))
{
arma::mat trainingData = std::move(CLI::GetParam<arma::mat>("training"));
// Cross-validation here.
size_t folds = CLI::GetParam<int>("folds");
if (folds == 0)
{
folds = trainingData.n_cols;
Log::Info << "Performing leave-one-out cross validation." << endl;
}
else
{
Log::Info << "Performing " << folds << "-fold cross validation." << endl;
}
trainingData = std::move(CLI::GetParam<arma::mat>("training"));
const bool regularization = false;
// const bool regularization = CLI::HasParam("volume_regularization");
const int maxLeafSize = CLI::GetParam<int>("max_leaf_size");
const int minLeafSize = CLI::GetParam<int>("min_leaf_size");
const bool skipPruning = CLI::HasParam("skip_pruning");
size_t folds = CLI::GetParam<int>("folds");
if (folds == 0)
folds = trainingData.n_cols;
// Obtain the optimal tree.
Timer::Start("det_training");
tree = Trainer<arma::mat, int>(trainingData, folds, regularization,
maxLeafSize, minLeafSize, "");
maxLeafSize, minLeafSize,
skipPruning);
Timer::Stop("det_training");
// Compute training set estimates, if desired.
@@ -166,25 +197,112 @@ void mlpackMain()
// the given file.
if (CLI::HasParam("test"))
{
arma::mat testData = std::move(CLI::GetParam<arma::mat>("test"));
// Compute test set densities.
Timer::Start("det_test_set_estimation");
arma::rowvec testDensities(testData.n_cols);
for (size_t i = 0; i < testData.n_cols; i++)
testDensities[i] = tree->ComputeValue(testData.unsafe_col(i));
Timer::Stop("det_test_set_estimation");
testData = std::move(CLI::GetParam<arma::mat>("test"));
if (CLI::HasParam("test_set_estimates"))
{
// Compute test set densities.
Timer::Start("det_test_set_estimation");
arma::rowvec testDensities(testData.n_cols);
for (size_t i = 0; i < testData.n_cols; i++)
testDensities[i] = tree->ComputeValue(testData.unsafe_col(i));
Timer::Stop("det_test_set_estimation");
CLI::GetParam<arma::mat>("test_set_estimates") = std::move(testDensities);
}
// Print variable importance.
if (CLI::HasParam("vi"))
{
arma::vec importances;
tree->ComputeVariableImportance(importances);
CLI::GetParam<arma::mat>("vi") = importances.t();
}
}
// Print variable importance.
if (CLI::HasParam("vi"))
if (CLI::HasParam("tag_file"))
{
arma::vec importances;
tree->ComputeVariableImportance(importances);
CLI::GetParam<arma::mat>("vi") = std::move(importances.t());
const arma::mat& estimationData =
CLI::HasParam("test") ? testData : trainingData;
const string tagFile = CLI::GetParam<string>("tag_file");
std::ofstream ofs;
ofs.open(tagFile, std::ofstream::out);
arma::Row<size_t> counters;
Timer::Start("det_test_set_tagging");
if (!ofs.is_open())
{
Log::Warn << "Unable to open file '" << tagFile
<< "' to save tag membership info."
<< std::endl;
}
else if (CLI::HasParam("path_format"))
{
const bool reqCounters = CLI::HasParam("tag_counters_file");
const string pathFormat = CLI::GetParam<string>("path_format");
PathCacher::PathFormat theFormat;
if (pathFormat == "lr" || pathFormat == "LR")
theFormat = PathCacher::FormatLR;
else if (pathFormat == "lr-id" || pathFormat == "LR-ID")
theFormat = PathCacher::FormatLR_ID;
else if (pathFormat == "id-lr" || pathFormat == "ID-LR")
theFormat = PathCacher::FormatID_LR;
else
{
Log::Warn << "Unknown path format specified: '" << pathFormat
<< "'. Valid are: lr | lr-id | id-lr. Defaults to 'lr'." << endl;
theFormat = PathCacher::FormatLR;
}
PathCacher path(theFormat, tree);
counters.zeros(path.NumNodes());
for (size_t i = 0; i < estimationData.n_cols; i++)
{
int tag = tree->FindBucket(estimationData.unsafe_col(i));
ofs << tag << " " << path.PathFor(tag) << std::endl;
for (; tag >= 0 && reqCounters; tag = path.ParentOf(tag))
counters(tag) += 1;
}
ofs.close();
if (reqCounters)
{
ofs.open(CLI::GetParam<string>("tag_counters_file"),
std::ofstream::out);
for (size_t j = 0; j < counters.n_elem; ++j)
ofs << j << " "
<< counters(j) << " "
<< path.PathFor(j) << endl;
ofs.close();
}
}
else
{
int numLeaves = tree->TagTree();
counters.zeros(numLeaves);
for (size_t i = 0; i < estimationData.n_cols; i++)
{
const int tag = tree->FindBucket(estimationData.unsafe_col(i));
ofs << tag << std::endl;
counters(tag) += 1;
}
if (CLI::HasParam("tag_counters_file"))
data::Save(CLI::GetParam<string>("tag_counters_file"), counters);
}
Timer::Stop("det_test_set_tagging");
ofs.close();
}
// Save the model, if desired.
src/mlpack/methods/det/dt_utils.hpp
+74 -2
View File
@@ -35,7 +35,7 @@ void PrintLeafMembership(DTree<MatType, TagType>* dtree,
const MatType& data,
const arma::Mat<size_t>& labels,
const size_t numClasses,
const std::string leafClassMembershipFile = "");
const std::string& leafClassMembershipFile = "");
/**
* Print the variable importance of each dimension of a density estimation tree.
@@ -67,7 +67,79 @@ DTree<MatType, TagType>* Trainer(MatType& dataset,
const bool useVolumeReg = false,
const size_t maxLeafSize = 10,
const size_t minLeafSize = 5,
const std::string unprunedTreeOutput = "");
const std::string unprunedTreeOutput = "",
const bool skipPruning = false);
/**
* This class is responsible for caching the path to each node of the tree. Its
* instance is provided to EnumerateTree() utility ONCE and it caches the paths
* to all the leafs and then easily (and quickly) retrieves these paths for each
* test entry.
*/
class PathCacher
{
public:
/**
* Possible formats to use for output.
*/
enum PathFormat
{
//! Print only whether we went left or right.
FormatLR,
//! Print the direction, then the tag of the node.
FormatLR_ID,
//! Print the tag of the node, then the direction.
FormatID_LR
};
/**
* Construct a PathCacher object on the given tree with the given format.
*
* @param fmt Format to use for output.
* @param tree Tree to cache paths in.
*/
template<typename MatType>
PathCacher(PathFormat fmt, DTree<MatType, int>* tree);
/**
* Enter a given node.
*/
template<typename MatType>
void Enter(const DTree<MatType, int>* node,
const DTree<MatType, int>* parent);
/**
* Leave the given node.
*/
template<typename MatType>
void Leave(const DTree<MatType, int>* node,
const DTree<MatType, int>* parent);
/**
* Return the constructed path for a given tag.
*/
const std::string& PathFor(int tag) const;
/**
* Get the parent tag of a given tag.
*/
int ParentOf(int tag) const;
/**
* Get the number of nodes in the path cache.
*/
size_t NumNodes() const { return pathCache.size(); }
protected:
typedef std::list<std::pair<bool, int>> PathType;
typedef std::vector<std::pair<int, std::string>> PathCacheType;
PathType path;
PathFormat format;
PathCacheType pathCache;
std::string BuildString();
};
} // namespace det
} // namespace mlpack
src/mlpack/methods/det/dt_utils_impl.hpp
+107 -46
View File
@@ -14,19 +14,20 @@
#define MLPACK_METHODS_DET_DT_UTILS_IMPL_HPP
#include "dt_utils.hpp"
#include <mlpack/core/tree/enumerate_tree.hpp>
namespace mlpack {
namespace det {
template <typename MatType, typename TagType>
void PrintLeafMembership(DTree<MatType, TagType>* dtree,
template <typename MatType>
void PrintLeafMembership(DTree<MatType, int>* dtree,
const MatType& data,
const arma::Mat<size_t>& labels,
const size_t numClasses,
const std::string leafClassMembershipFile)
const std::string& leafClassMembershipFile)
{
// Tag the leaves with numbers.
TagType numLeaves = dtree->TagTree();
int numLeaves = dtree->TagTree();
arma::Mat<size_t> table(numLeaves, (numClasses + 1));
table.zeros();
@@ -34,7 +35,7 @@ void PrintLeafMembership(DTree<MatType, TagType>* dtree,
for (size_t i = 0; i < data.n_cols; i++)
{
const typename MatType::vec_type testPoint = data.unsafe_col(i);
const TagType leafTag = dtree->FindBucket(testPoint);
const int leafTag = dtree->FindBucket(testPoint);
const size_t label = labels[i];
table(leafTag, label) += 1;
}
@@ -111,11 +112,12 @@ DTree<MatType, TagType>* Trainer(MatType& dataset,
const bool useVolumeReg,
const size_t maxLeafSize,
const size_t minLeafSize,
const std::string unprunedTreeOutput)
const bool skipPruning)
{
// Initialize the tree.
DTree<MatType, TagType> dtree(dataset);
DTree<MatType, TagType>* dtree = new DTree<MatType, TagType>(dataset);
Timer::Start("tree_growing");
// Prepare to grow the tree...
arma::Col<size_t> oldFromNew(dataset.n_cols);
for (size_t i = 0; i < oldFromNew.n_elem; i++)
@@ -126,60 +128,52 @@ DTree<MatType, TagType>* Trainer(MatType& dataset,
// Growing the tree
double oldAlpha = 0.0;
double alpha = dtree.Grow(newDataset, oldFromNew, useVolumeReg, maxLeafSize,
double alpha = dtree->Grow(newDataset, oldFromNew, useVolumeReg, maxLeafSize,
minLeafSize);
Log::Info << dtree.SubtreeLeaves() << " leaf nodes in the tree using full "
Timer::Stop("tree_growing");
Log::Info << dtree->SubtreeLeaves() << " leaf nodes in the tree using full "
<< "dataset; minimum alpha: " << alpha << "." << std::endl;
// Compute densities for the training points in the full tree, if we were
// asked for this.
if (unprunedTreeOutput != "")
{
std::ofstream outfile(unprunedTreeOutput.c_str());
if (outfile.good())
{
for (size_t i = 0; i < dataset.n_cols; ++i)
{
arma::vec testPoint = dataset.unsafe_col(i);
outfile << dtree.ComputeValue(testPoint) << std::endl;
}
}
else
{
Log::Warn << "Can't open '" << unprunedTreeOutput << "' to write computed"
<< " densities to." << std::endl;
}
if (skipPruning)
return dtree;
outfile.close();
}
if (folds == dataset.n_cols)
Log::Info << "Performing leave-one-out cross validation." << std::endl;
else
Log::Info << "Performing " << folds << "-fold cross validation." <<
std::endl;
Timer::Start("pruning_sequence");
// Sequentially prune and save the alpha values and the values of c_t^2 * r_t.
std::vector<std::pair<double, double> > prunedSequence;
while (dtree.SubtreeLeaves() > 1)
while (dtree->SubtreeLeaves() > 1)
{
std::pair<double, double> treeSeq(oldAlpha,
dtree.SubtreeLeavesLogNegError());
dtree->SubtreeLeavesLogNegError());
prunedSequence.push_back(treeSeq);
oldAlpha = alpha;
alpha = dtree.PruneAndUpdate(oldAlpha, dataset.n_cols, useVolumeReg);
alpha = dtree->PruneAndUpdate(oldAlpha, dataset.n_cols, useVolumeReg);
// Some sanity checks. It seems that on some datasets, the error does not
// increase as the tree is pruned but instead stays the same---hence the
// "<=" in the final assert.
Log::Assert((alpha < std::numeric_limits<double>::max())
|| (dtree.SubtreeLeaves() == 1));
|| (dtree->SubtreeLeaves() == 1));
Log::Assert(alpha > oldAlpha);
Log::Assert(dtree.SubtreeLeavesLogNegError() <= treeSeq.second);
Log::Assert(dtree->SubtreeLeavesLogNegError() <= treeSeq.second);
}
std::pair<double, double> treeSeq(oldAlpha, dtree.SubtreeLeavesLogNegError());
std::pair<double, double> treeSeq(oldAlpha,
dtree->SubtreeLeavesLogNegError());
prunedSequence.push_back(treeSeq);
Timer::Stop("pruning_sequence");
Log::Info << prunedSequence.size() << " trees in the sequence; maximum alpha:"
<< " " << oldAlpha << "." << std::endl;
MatType cvData(dataset);
const MatType cvData(dataset);
const size_t testSize = dataset.n_cols / folds;
arma::vec regularizationConstants(prunedSequence.size());
@@ -191,7 +185,7 @@ DTree<MatType, TagType>* Trainer(MatType& dataset,
// implementation. omp_size_t is the appropriate type according to the
// platform.
#pragma omp parallel for default(none) \
shared(cvData, prunedSequence, regularizationConstants)
shared(prunedSequence, regularizationConstants)
for (omp_size_t fold = 0; fold < (omp_size_t) folds; fold++)
{
// Break up data into train and test sets.
@@ -289,8 +283,9 @@ DTree<MatType, TagType>* Trainer(MatType& dataset,
Log::Info << "Optimal alpha: " << optimalAlpha << "." << std::endl;
// Initialize the tree.
DTree<MatType, TagType>* dtreeOpt = new DTree<MatType, TagType>(dataset);
// Re-Initialize the tree.
delete dtree;
dtree = new DTree<MatType, TagType>(dataset);
// Getting ready to grow the tree...
for (size_t i = 0; i < oldFromNew.n_elem; i++)
@@ -301,31 +296,97 @@ DTree<MatType, TagType>* Trainer(MatType& dataset,
// Grow the tree.
oldAlpha = -DBL_MAX;
alpha = dtreeOpt->Grow(newDataset,
alpha = dtree->Grow(newDataset,
oldFromNew,
useVolumeReg,
maxLeafSize,
minLeafSize);
// Prune with optimal alpha.
while ((oldAlpha < optimalAlpha) && (dtreeOpt->SubtreeLeaves() > 1))
while ((oldAlpha < optimalAlpha) && (dtree->SubtreeLeaves() > 1))
{
oldAlpha = alpha;
alpha = dtreeOpt->PruneAndUpdate(oldAlpha, newDataset.n_cols, useVolumeReg);
alpha = dtree->PruneAndUpdate(oldAlpha, newDataset.n_cols, useVolumeReg);
// Some sanity checks.
Log::Assert((alpha < std::numeric_limits<double>::max()) ||
(dtreeOpt->SubtreeLeaves() == 1));
(dtree->SubtreeLeaves() == 1));
Log::Assert(alpha > oldAlpha);
}
Log::Info << dtreeOpt->SubtreeLeaves() << " leaf nodes in the optimally "
Log::Info << dtree->SubtreeLeaves() << " leaf nodes in the optimally "
<< "pruned tree; optimal alpha: " << oldAlpha << "." << std::endl;
return dtreeOpt;
return dtree;
}
template<typename MatType>
PathCacher::PathCacher(PathCacher::PathFormat fmt, DTree<MatType, int>* dtree) :
format(fmt)
{
// Here we use TagTree()'s output to determine the
// number of _nodes_ in the tree.
pathCache.resize(dtree->TagTree(0, true));
pathCache[0] = PathCacheType::value_type(-1, "");
tree::EnumerateTree(dtree, *this);
}
template<typename MatType>
void PathCacher::Enter(const DTree<MatType, int>* node,
const DTree<MatType, int>* parent)
{
if (parent == nullptr)
return;
int tag = node->BucketTag();
path.push_back(PathType::value_type(parent->Left() == node, tag));
pathCache[tag] = PathCacheType::value_type(parent->BucketTag(),
(node->SubtreeLeaves() > 1) ?
"" : BuildString());
}
template<typename MatType>
void PathCacher::Leave(const DTree<MatType, int>* /* node */,
const DTree<MatType, int>* parent)
{
if (parent != nullptr)
path.pop_back();
}
std::string PathCacher::BuildString()
{
std::string str("");
for (PathType::iterator it = path.begin(); it != path.end(); it++)
{
switch (format)
{
case FormatLR:
str += it->first ? "L" : "R";
break;
case FormatLR_ID:
str += (it->first ? "L" : "R") + std::to_string(it->second);
break;
case FormatID_LR:
str += std::to_string(it->second) + (it->first ? "L" : "R");
break;
}
}
return str;
}
int PathCacher::ParentOf(int tag) const
{
return pathCache[tag].first;
}
const std::string& PathCacher::PathFor(int tag) const
{
return pathCache[tag].second;
}
} // namespace det
} // namespace mlpack
#endif
#endif // MLPACK_METHODS_DET_DT_UTILS_IMPL_HPP
src/mlpack/methods/det/dtree.hpp
+28 -9
View File
@@ -46,12 +46,12 @@ template<typename MatType = arma::mat,
class DTree
{
public:
/**
* The actual, underlying type we're working with.
*/
typedef typename MatType::elem_type ElemType;
typedef typename MatType::vec_type VecType;
typedef typename arma::Col<ElemType> StatType;
//! The actual, underlying type we're working with.
typedef typename MatType::elem_type ElemType;
//! The type of vector we are using.
typedef typename MatType::vec_type VecType;
//! The statistic type we are holding.
typedef typename arma::Col<ElemType> StatType;
/**
* Create an empty density estimation tree.
@@ -184,11 +184,14 @@ class DTree
/**
* Index the buckets for possible usage later; this results in every leaf in
* the tree having a specific tag (accessible with BucketTag()). This
* function calls itself recursively.
* function calls itself recursively. The tag is incremented with
* `operator++()`, so any `TagType` overriding it will do.
*
* @param tag Tag for the next leaf; leave at 0 for the initial call.
* @param everyNodde Whether to increment on every node, not just leaves.
*/
TagType TagTree(const TagType& tag = 0);
TagType TagTree(const TagType& tag = 0, bool everyNode = false);
/**
* Return the tag of the leaf containing the query. This is useful for
@@ -198,6 +201,7 @@ class DTree
*/
TagType FindBucket(const VecType& query) const;
/**
* Compute the variable importance of each dimension in the learned tree.
*
@@ -301,7 +305,19 @@ class DTree
//! Return the upper part of the alpha sum.
double AlphaUpper() const { return alphaUpper; }
//! Return the current bucket's ID, if leaf, or -1 otherwise
TagType BucketTag() const { return subtreeLeaves == 1 ? bucketTag : -1; }
TagType BucketTag() const { return bucketTag; }
//! Return the number of children in this node.
size_t NumChildren() const { return !left ? 0 : 2; }
/**
* Return the specified child (0 will be left, 1 will be right). If the index
* is greater than 1, this will return the right child.
*
* @param child Index of child to return.
*/
DTree& Child(const size_t child) const { return !child ? *left : *right; }
DTree*& ChildPtr(const size_t child) { return (!child) ? left : right; }
//! Return the maximum values.
const StatType& MaxVals() const { return maxVals; }
@@ -335,6 +351,9 @@ class DTree
const size_t splitDim,
const ElemType splitValue,
arma::Col<size_t>& oldFromNew) const;
void FillMinMax(const StatType& mins,
const StatType& maxs);
};
} // namespace det
src/mlpack/methods/det/dtree_impl.hpp
+206 -153
View File
@@ -20,134 +20,136 @@ using namespace det;
namespace details
{
/**
* This one sorts and scand the given per-dimension extract and puts all splits
* in a vector, that can easily be iterated afterwards. General implementation.
*/
template <typename ElemType, typename MatType>
void ExtractSplits(std::vector<std::pair<ElemType, size_t>>& splitVec,
const MatType& data,
size_t dim,
const size_t start,
const size_t end,
const size_t minLeafSize)
{
static_assert(
std::is_same<typename MatType::elem_type, ElemType>::value == true,
"The ElemType does not correspond to the matrix's element type.");
typedef std::pair<ElemType, size_t> SplitItem;
const typename MatType::row_type dimVec =
/**
* This one sorts and scand the given per-dimension extract and puts all splits
* in a vector, that can easily be iterated afterwards. General implementation.
*/
template<typename ElemType, typename MatType>
void ExtractSplits(std::vector<std::pair<ElemType, size_t>>& splitVec,
const MatType& data,
size_t dim,
const size_t start,
const size_t end,
const size_t minLeafSize)
{
static_assert(
std::is_same<typename MatType::elem_type, ElemType>::value == true,
"The ElemType does not correspond to the matrix's element type.");
typedef std::pair<ElemType, size_t> SplitItem;
const typename MatType::row_type dimVec =
arma::sort(data(dim, arma::span(start, end - 1)));
// Ensure the minimum leaf size on both sides. We need to figure out why
// there are spikes if this minLeafSize is enforced here...
for (size_t i = minLeafSize - 1; i < dimVec.n_elem - minLeafSize; ++i)
// Ensure the minimum leaf size on both sides. We need to figure out why there
// are spikes if this minLeafSize is enforced here...
for (size_t i = minLeafSize - 1; i < dimVec.n_elem - minLeafSize; ++i)
{
// This makes sense for real continuous data. This kinda corrupts the data
// and estimation if the data is ordinal. Potentially we can fix that by
// taking into account ordinality later in the min/max update, but then we
// can end-up with a zero-volumed dimension. No good.
const ElemType split = (dimVec[i] + dimVec[i + 1]) / 2.0;
// Check if we can split here (two points are different)
if (split != dimVec[i])
splitVec.push_back(SplitItem(split, i + 1));
}
}
// Now the custom arma::Mat implementation.
template<typename ElemType>
void ExtractSplits(std::vector<std::pair<ElemType, size_t>>& splitVec,
const arma::Mat<ElemType>& data,
size_t dim,
const size_t start,
const size_t end,
const size_t minLeafSize)
{
typedef std::pair<ElemType, size_t> SplitItem;
arma::rowvec dimVec = data(dim, arma::span(start, end - 1));
// We sort these, in-place (it's a copy of the data, anyways).
std::sort(dimVec.begin(), dimVec.end());
for (size_t i = minLeafSize - 1; i < dimVec.n_elem - minLeafSize; ++i)
{
// This makes sense for real continuous data. This kinda corrupts the data
// and estimation if the data is ordinal. Potentially we can fix that by
// taking into account ordinality later in the min/max update, but then we
// can end-up with a zero-volumed dimension. No good.
const ElemType split = (dimVec[i] + dimVec[i + 1]) / 2.0;
if (split != dimVec[i])
splitVec.push_back(SplitItem(split, i + 1));
}
}
// This the custom, sparse optimized implementation of the same routine.
template<typename ElemType>
void ExtractSplits(std::vector<std::pair<ElemType, size_t>>& splitVec,
const arma::SpMat<ElemType>& data,
size_t dim,
const size_t start,
const size_t end,
const size_t minLeafSize)
{
// It's common sense, but we also use it in a check later.
Log::Assert(minLeafSize > 0);
typedef std::pair<ElemType, size_t> SplitItem;
const size_t n_elem = end - start;
// Construct a vector of values.
const arma::SpRow<ElemType> row = data(dim, arma::span(start, end - 1));
std::vector<ElemType> valsVec(row.begin(), row.end());
// ... and sort it!
std::sort(valsVec.begin(), valsVec.end());
// Now iterate over the values, taking account for the over-the-zeroes jump
// and construct the splits vector.
const size_t zeroes = n_elem - valsVec.size();
ElemType lastVal = -std::numeric_limits<ElemType>::max();
size_t padding = 0;
for (size_t i = 0; i < valsVec.size(); ++i)
{
const ElemType newVal = valsVec[i];
if (lastVal < ElemType(0) && newVal > ElemType(0) && zeroes > 0)
{
// This makes sense for real continuous data. This kinda corrupts the
// data and estimation if the data is ordinal. Potentially we can fix
// that by taking into account ordinality later in the min/max update,
// but then we can end-up with a zero-volumed dimension. No good.
const ElemType split = (dimVec[i] + dimVec[i + 1]) / 2.0;
Log::Assert(padding == 0); // We should arrive here once!
// The minLeafSize > 0 also guarantees we're not entering right at the
// start.
if (i >= minLeafSize && i <= n_elem - minLeafSize)
splitVec.push_back(SplitItem(lastVal / 2.0, i));
padding = zeroes;
lastVal = ElemType(0);
}
// This is the normal case.
if (i + padding >= minLeafSize && i + padding <= n_elem - minLeafSize)
{
// This makes sense for real continuous data. This kinda corrupts the
// data and estimation if the data is ordinal. Potentially we can fix that
// by taking into account ordinality later in the min/max update, but then
// we can end-up with a zero-volumed dimension. No good.
const ElemType split = (lastVal + newVal) / 2.0;
// Check if we can split here (two points are different)
if (split != dimVec[i])
splitVec.push_back(SplitItem(split, i + 1));
if (split != newVal)
splitVec.push_back(SplitItem(split, i + padding));
}
lastVal = newVal;
}
}
// Now the custom arma::Mat implementation
template <typename ElemType>
void ExtractSplits(std::vector<std::pair<ElemType, size_t>>& splitVec,
const arma::Mat<ElemType>& data,
size_t dim,
const size_t start,
const size_t end,
const size_t minLeafSize)
{
typedef std::pair<ElemType, size_t> SplitItem;
arma::vec dimVec = data(dim, arma::span(start, end - 1)).t();
} // namespace details
// We sort these, in-place (it's a copy of the data, anyways).
std::sort(dimVec.begin(), dimVec.end());
for (size_t i = minLeafSize - 1; i < dimVec.n_elem - minLeafSize; ++i)
{
// This makes sense for real continuous data. This kinda corrupts the
// data and estimation if the data is ordinal. Potentially we can fix
// that by taking into account ordinality later in the min/max update,
// but then we can end-up with a zero-volumed dimension. No good.
const ElemType split = (dimVec[i] + dimVec[i + 1]) / 2.0;
if (split != dimVec[i])
splitVec.push_back(SplitItem(split, i + 1));
}
}
// This the custom, sparse optimized implementation of the same routine.
template <typename ElemType>
void ExtractSplits(std::vector<std::pair<ElemType, size_t>>& splitVec,
const arma::SpMat<ElemType>& data,
size_t dim,
const size_t start,
const size_t end,
const size_t minLeafSize)
{
// It's common sense, but we also use it in a check later.
Log::Assert(minLeafSize > 0);
typedef std::pair<ElemType, size_t> SplitItem;
const size_t n_elem = end - start;
// Construct a vector of values.
const arma::SpRow<ElemType> row = data(dim, arma::span(start, end - 1));
std::vector<ElemType> valsVec(row.begin(), row.end());
// ... and sort it!
std::sort(valsVec.begin(), valsVec.end());
// Now iterate over the values, taking account for the over-the-zeroes
// jump and construct the splits vector.
const size_t zeroes = n_elem - valsVec.size();
ElemType lastVal = -std::numeric_limits<ElemType>::max();
size_t padding = 0;
for (size_t i = 0; i < valsVec.size(); ++i)
{
const ElemType newVal = valsVec[i];
if (lastVal < ElemType(0) && newVal > ElemType(0) && zeroes > 0)
{
Log::Assert(padding == 0); // We should arrive here once!
// The minLeafSize > 0 also guarantees we're not entering right at the
// start.
if (i >= minLeafSize && i <= n_elem - minLeafSize)
splitVec.push_back(SplitItem(lastVal / 2.0, i));
padding = zeroes;
lastVal = ElemType(0);
}
// the normal case
if (i + padding >= minLeafSize && i + padding <= n_elem - minLeafSize)
{
// This makes sense for real continuous data. This kinda corrupts the
// data and estimation if the data is ordinal. Potentially we can fix
// that by taking into account ordinality later in the min/max update,
// but then we can end-up with a zero-volumed dimension. No good.
const ElemType split = (lastVal + newVal) / 2.0;
// Check if we can split here (two points are different)
if (split != newVal)
splitVec.push_back(SplitItem(split, i + padding));
}
lastVal = newVal;
}
}
}; // namespace details
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
DTree<MatType, TagType>::DTree() :
start(0),
end(0),
@@ -165,7 +167,7 @@ DTree<MatType, TagType>::DTree() :
right(NULL)
{ /* Nothing to do. */ }
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
DTree<MatType, TagType>::DTree(const DTree& obj) :
start(obj.start),
end(obj.end),
@@ -187,10 +189,13 @@ DTree<MatType, TagType>::DTree(const DTree& obj) :
/* Nothing to do. */
}
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
DTree<MatType, TagType>& DTree<MatType, TagType>::operator=(
const DTree<MatType, TagType>& obj)
{
if (this == &obj)
return *this;
// Copy the values from the other tree.
start = obj.start;
end = obj.end;
@@ -213,12 +218,12 @@ DTree<MatType, TagType>& DTree<MatType, TagType>::operator=(
// Copy the children.
left = ((obj.left == NULL) ? NULL : new DTree(*obj.left));
left = ((obj.right == NULL) ? NULL : new DTree(*obj.right));
right = ((obj.right == NULL) ? NULL : new DTree(*obj.right));
return *this;
}
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
DTree<MatType, TagType>::DTree(DTree&& obj):
start(obj.start),
end(obj.end),
@@ -254,10 +259,13 @@ DTree<MatType, TagType>::DTree(DTree&& obj):
obj.right = NULL;
}
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
DTree<MatType, TagType>& DTree<MatType, TagType>::operator=(
DTree<MatType, TagType>&& obj)
{
if (this == &obj)
return *this;
// Move the values from the other tree.
start = obj.start;
end = obj.end;
@@ -302,8 +310,8 @@ DTree<MatType, TagType>& DTree<MatType, TagType>::operator=(
}
// Root node initializers
template <typename MatType, typename TagType>
// Root node initializers.
template<typename MatType, typename TagType>
DTree<MatType, TagType>::DTree(const StatType& maxVals,
const StatType& minVals,
const size_t totalPoints) :
@@ -325,7 +333,7 @@ DTree<MatType, TagType>::DTree(const StatType& maxVals,
right(NULL)
{ /* Nothing to do. */ }
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
DTree<MatType, TagType>::DTree(MatType & data) :
start(0),
end(data.n_cols),
@@ -346,8 +354,8 @@ DTree<MatType, TagType>::DTree(MatType & data) :
logNegError = LogNegativeError(data.n_cols);
}
// Non-root node initializers
template <typename MatType, typename TagType>
// Non-root node initializers.
template<typename MatType, typename TagType>
DTree<MatType, TagType>::DTree(const StatType& maxVals,
const StatType& minVals,
const size_t start,
@@ -371,7 +379,7 @@ DTree<MatType, TagType>::DTree(const StatType& maxVals,
right(NULL)
{ /* Nothing to do. */ }
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
DTree<MatType, TagType>::DTree(const StatType& maxVals,
const StatType& minVals,
const size_t totalPoints,
@@ -395,7 +403,7 @@ DTree<MatType, TagType>::DTree(const StatType& maxVals,
right(NULL)
{ /* Nothing to do. */ }
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
DTree<MatType, TagType>::~DTree()
{
delete left;
@@ -404,7 +412,7 @@ DTree<MatType, TagType>::~DTree()
// This function computes the log-l2-negative-error of a given node from the
// formula R(t) = log(|t|^2 / (N^2 V_t)).
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
double DTree<MatType, TagType>::LogNegativeError(const size_t totalPoints) const
{
// log(-|t|^2 / (N^2 V_t)) = log(-1) + 2 log(|t|) - 2 log(N) - log(V_t).
@@ -425,7 +433,7 @@ double DTree<MatType, TagType>::LogNegativeError(const size_t totalPoints) const
// This function finds the best split with respect to the L2-error, by trying
// all possible splits. The dataset is the full data set but the start and
// end are used to obtain the point in this node.
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
bool DTree<MatType, TagType>::FindSplit(const MatType& data,
size_t& splitDim,
ElemType& splitValue,
@@ -433,7 +441,7 @@ bool DTree<MatType, TagType>::FindSplit(const MatType& data,
double& rightError,
const size_t minLeafSize) const
{
typedef std::pair<ElemType, size_t> SplitItem;
typedef std::pair<ElemType, size_t> SplitItem;
// Ensure the dimensionality of the data is the same as the dimensionality of
// the bounding rectangle.
@@ -543,7 +551,7 @@ bool DTree<MatType, TagType>::FindSplit(const MatType& data,
return splitFound;
}
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
size_t DTree<MatType, TagType>::SplitData(MatType& data,
const size_t splitDim,
const ElemType splitValue,
@@ -577,8 +585,8 @@ size_t DTree<MatType, TagType>::SplitData(MatType& data,
return left;
}
// Greedily expand the tree
template <typename MatType, typename TagType>
// Greedily expand the tree.
template<typename MatType, typename TagType>
double DTree<MatType, TagType>::Grow(MatType& data,
arma::Col<size_t>& oldFromNew,
const bool useVolReg,
@@ -727,7 +735,7 @@ double DTree<MatType, TagType>::Grow(MatType& data,
}
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
double DTree<MatType, TagType>::PruneAndUpdate(const double oldAlpha,
const size_t points,
const bool useVolReg)
@@ -842,7 +850,7 @@ double DTree<MatType, TagType>::PruneAndUpdate(const double oldAlpha,
//
// Future improvement: Open up the range with epsilons on both sides where
// epsilon depends on the density near the boundary.
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
bool DTree<MatType, TagType>::WithinRange(const VecType& query) const
{
for (size_t i = 0; i < query.n_elem; ++i)
@@ -853,7 +861,7 @@ bool DTree<MatType, TagType>::WithinRange(const VecType& query) const
}
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
double DTree<MatType, TagType>::ComputeValue(const VecType& query) const
{
Log::Assert(query.n_elem == maxVals.n_elem);
@@ -882,8 +890,8 @@ double DTree<MatType, TagType>::ComputeValue(const VecType& query) const
}
// Index the buckets for possible usage later.
template <typename MatType, typename TagType>
TagType DTree<MatType, TagType>::TagTree(const TagType& tag)
template<typename MatType, typename TagType>
TagType DTree<MatType, TagType>::TagTree(const TagType& tag, bool every)
{
if (subtreeLeaves == 1)
{
@@ -891,18 +899,33 @@ TagType DTree<MatType, TagType>::TagTree(const TagType& tag)
bucketTag = tag;
return (tag + 1);
}
else
TagType nextTag;
if (every)
{
return right->TagTree(left->TagTree(tag));
bucketTag = tag;
nextTag = (tag + 1);
}
else
nextTag = tag;
return right->TagTree(left->TagTree(nextTag, every), every);
}
template <typename MatType, typename TagType>
template<typename MatType, typename TagType>
TagType DTree<MatType, TagType>::FindBucket(const VecType& query) const
{
Log::Assert(query.n_elem == maxVals.n_elem);
if (subtreeLeaves == 1) // If we are a leaf...
if (root == 1) // If we are the root...
{
// Check if the query is within range.
if (!WithinRange(query))
return -1;
}
// If we are a leaf...
if (subtreeLeaves == 1)
{
return bucketTag;
}
@@ -915,9 +938,9 @@ TagType DTree<MatType, TagType>::FindBucket(const VecType& query) const
}
}
template <typename MatType, typename TagType>
void
DTree<MatType, TagType>::ComputeVariableImportance(arma::vec& importances) const
template<typename MatType, typename TagType>
void DTree<MatType, TagType>::ComputeVariableImportance(arma::vec& importances)
const
{
// Clear and set to right size.
importances.zeros(maxVals.n_elem);
@@ -944,8 +967,31 @@ DTree<MatType, TagType>::ComputeVariableImportance(arma::vec& importances) const
}
}
template <typename MatType, typename TagType>
template <typename Archive>
template<typename MatType, typename TagType>
void DTree<MatType, TagType>::FillMinMax(const StatType& mins,
const StatType& maxs)
{
if (!root)
{
minVals = mins;
maxVals = maxs;
}
if (left && right)
{
StatType maxValsL(maxs);
StatType maxValsR(maxs);
StatType minValsL(mins);
StatType minValsR(mins);
maxValsL[splitDim] = minValsR[splitDim] = splitValue;
left->FillMinMax(minValsL, maxValsL);
right->FillMinMax(minValsR, maxValsR);
}
}
template<typename MatType, typename TagType>
template<typename Archive>
void DTree<MatType, TagType>::Serialize(Archive& ar,
const unsigned int /* version */)
{
@@ -953,8 +999,6 @@ void DTree<MatType, TagType>::Serialize(Archive& ar,
ar & CreateNVP(start, "start");
ar & CreateNVP(end, "end");
ar & CreateNVP(maxVals, "maxVals");
ar & CreateNVP(minVals, "minVals");
ar & CreateNVP(splitDim, "splitDim");
ar & CreateNVP(splitValue, "splitValue");
ar & CreateNVP(logNegError, "logNegError");
@@ -976,5 +1020,14 @@ void DTree<MatType, TagType>::Serialize(Archive& ar,
ar & CreateNVP(left, "left");
ar & CreateNVP(right, "right");
}
if (root)
{
ar & CreateNVP(maxVals, "maxVals");
ar & CreateNVP(minVals, "minVals");
// This is added in order to reduce (dramatically!) the model file size.
if (Archive::is_loading::value && left && right)
FillMinMax(minVals, maxVals);
}
}
src/mlpack/methods/gmm/diagonal_constraint.hpp
+1 -2
View File
@@ -27,8 +27,7 @@ class DiagonalConstraint
static void ApplyConstraint(arma::mat& covariance)
{
// Save the diagonal only.
arma::vec diagonal = covariance.diag();
covariance = arma::diagmat(diagonal);
covariance = arma::diagmat(arma::clamp(covariance.diag(), 1e-10, DBL_MAX));
}
//! Serialize the constraint (which holds nothing, so, nothing to do).
src/mlpack/methods/gmm/em_fit.hpp
+20
View File
@@ -162,6 +162,26 @@ class EMFit
dists,
const arma::vec& weights) const;
// Armadillo uses uword internally as an OpenMP index type, which crashes
// Visual Studio.
#ifndef _WIN32
/**
* Use the Armadillo gmm_diag clusterer to train a GMM with diagonal
* covariance. If InitialClusteringType == kmeans::KMeans<>, this will use
* Armadillo's initialization also.
*
* @param observations Data to train on.
* @param dists Distributions to store model in.
* @param weights Prior weights.
* @param useInitialModel If true, the existing model will be used.
*/
void ArmadilloGMMWrapper(
const arma::mat& observations,
std::vector<distribution::GaussianDistribution>& dists,
arma::vec& weights,
const bool useInitialModel);
#endif
//! Maximum iterations of EM algorithm.
size_t maxIterations;
//! Tolerance for convergence of EM.
src/mlpack/methods/gmm/em_fit_impl.hpp
+74
View File
@@ -15,6 +15,7 @@
// In case it hasn't been included yet.
#include "em_fit.hpp"
#include "diagonal_constraint.hpp"
namespace mlpack {
namespace gmm {
@@ -39,6 +40,17 @@ void EMFit<InitialClusteringType, CovarianceConstraintPolicy>::Estimate(
arma::vec& weights,
const bool useInitialModel)
{
// Shortcut: if the user is using the DiagonalConstraint, then we will call
// out to Armadillo. But Armadillo uses uword internally as an OpenMP index
// type, which crashes Visual Studio, so don't do this on Windows.
#ifndef _WIN32
if (std::is_same<CovarianceConstraintPolicy, DiagonalConstraint>::value)
{
ArmadilloGMMWrapper(observations, dists, weights, useInitialModel);
return;
}
#endif
// Only perform initial clustering if the user wanted it.
if (!useInitialModel)
InitialClustering(observations, dists, weights);
@@ -316,6 +328,68 @@ void EMFit<InitialClusteringType, CovarianceConstraintPolicy>::Serialize(
ar & CreateNVP(constraint, "constraint");
}
// Armadillo uses uword internally as an OpenMP index type, which crashes Visual
// Studio.
#ifndef _WIN32
template<typename InitialClusteringType, typename CovarianceConstraintPolicy>
void EMFit<InitialClusteringType, CovarianceConstraintPolicy>::
ArmadilloGMMWrapper(const arma::mat& observations,
std::vector<distribution::GaussianDistribution>& dists,
arma::vec& weights,
const bool useInitialModel)
{
arma::gmm_diag g;
// Warn the user that tolerance isn't used for convergence here if they've
// specified a non-default value.
if (tolerance != EMFit().Tolerance())
Log::Warn << "GMM::Train(): tolerance ignored when training GMMs with "
<< "DiagonalConstraint." << std::endl;
// If the initial clustering is the default k-means, we'll just use
// Armadillo's implementation. If mlpack ever changes k-means defaults to use
// something that is reliably quicker than the Lloyd iteration k-means update,
// then this code maybe should be revisited.
if (!std::is_same<InitialClusteringType, mlpack::kmeans::KMeans<>>::value ||
useInitialModel)
{
// Use clusterer to get initial values.
if (!useInitialModel)
InitialClustering(observations, dists, weights);
// Assemble matrix of means.
arma::mat means(observations.n_rows, dists.size());
arma::mat covs(observations.n_rows, dists.size());
for (size_t i = 0; i < dists.size(); ++i)
{
means.col(i) = dists[i].Mean();
covs.col(i) = dists[i].Covariance().diag();
}
g.reset(observations.n_rows, dists.size());
g.set_params(std::move(means), std::move(covs), weights.t());
g.learn(observations, dists.size(), arma::eucl_dist, arma::keep_existing, 0,
maxIterations, 1e-10, false /* no printing */);
}
else
{
// Use Armadillo for the initial clustering. We'll try and match mlpack
// defaults.
g.learn(observations, dists.size(), arma::eucl_dist, arma::static_subset,
1000, maxIterations, 1e-10, false /* no printing */);
}
// Extract means, covariances, and weights.
weights = g.hefts.t();
for (size_t i = 0; i < dists.size(); ++i)
{
dists[i].Mean() = g.means.col(i);
dists[i].Covariance(std::move(arma::diagmat(g.dcovs.col(i))));
}
}
#endif
} // namespace gmm
} // namespace mlpack
src/mlpack/methods/gmm/gmm_train_main.cpp
+38 -5
View File
@@ -15,6 +15,7 @@
#include "gmm.hpp"
#include "no_constraint.hpp"
#include "diagonal_constraint.hpp"
#include <mlpack/methods/kmeans/refined_start.hpp>
@@ -52,6 +53,11 @@ PROGRAM_INFO("Gaussian Mixture Model (GMM) Training",
"Bradley-Fayyad refined start initialization will be used. This can often "
"lead to better clustering results."
"\n\n"
"The 'diagonal_covariance' flag will cause the learned covariances to be "
"diagonal matrices. This significantly simplifies the model itself and "
"causes training to be faster, but restricts the ability to fit more "
"complex GMMs."
"\n\n"
"If GMM training fails with an error indicating that a covariance matrix "
"could not be inverted, make sure that the " +
PRINT_PARAM_STRING("no_force_positive") + " parameter is not "
@@ -94,6 +100,8 @@ PARAM_FLAG("no_force_positive", "Do not force the covariance matrices to be "
"positive definite.", "P");
PARAM_INT_IN("max_iterations", "Maximum number of iterations of EM algorithm "
"(passing 0 will run until convergence).", "n", 250);
PARAM_FLAG("diagonal_covariance", "Force the covariance of the Gaussians to "
"be diagonal. This can accelerate training time significantly.", "d");
// Parameters for dataset modification.
PARAM_DOUBLE_IN("noise", "Variance of zero-mean Gaussian noise to add to data.",
@@ -128,6 +136,11 @@ void mlpackMain()
"be greater than or equal to 1." << std::endl;
}
if (CLI::HasParam("diagonal_covariance") &&
CLI::HasParam("no_force_positive"))
Log::Warn << "--no_force_positive ignored because --diagonal_covariance is "
<< "specified!" << endl;
if (!CLI::HasParam("output_model"))
Log::Warn << "--output_model_file is not specified, so no model will be "
<< "saved!" << endl;
@@ -163,6 +176,7 @@ void mlpackMain()
const size_t maxIterations = (size_t) CLI::GetParam<int>("max_iterations");
const double tolerance = CLI::GetParam<double>("tolerance");
const bool forcePositive = !CLI::HasParam("no_force_positive");
const bool diagonalCovariance = CLI::HasParam("diagonal_covariance");
// This gets a bit weird because we need different types depending on whether
// --refined_start is specified.
@@ -186,9 +200,18 @@ void mlpackMain()
KMeansType k(1000, metric::SquaredEuclideanDistance(),
RefinedStart(samplings, percentage));
// Depending on the value of 'forcePositive', we have to use different
// types.
if (forcePositive)
// Depending on the value of forcePositive and diagonalCovariance, we have
// to use different types.
if (diagonalCovariance)
{
// Compute the parameters of the model using the EM algorithm.
Timer::Start("em");
EMFit<KMeansType, DiagonalConstraint> em(maxIterations, tolerance, k);
likelihood = gmm.Train(dataPoints, CLI::GetParam<int>("trials"), false,
em);
Timer::Stop("em");
}
else if (forcePositive)
{
// Compute the parameters of the model using the EM algorithm.
Timer::Start("em");
@@ -209,8 +232,18 @@ void mlpackMain()
}
else
{
// Depending on the value of forcePositive, we have to use different types.
if (forcePositive)
// Depending on the value of forcePositive and diagonalCovariance, we have
// to use different types.
if (diagonalCovariance)
{
// Compute the parameters of the model using the EM algorithm.
Timer::Start("em");
EMFit<kmeans::KMeans<>, DiagonalConstraint> em(maxIterations, tolerance);
likelihood = gmm.Train(dataPoints, CLI::GetParam<int>("trials"), false,
em);
Timer::Stop("em");
}
else if (forcePositive)
{
// Compute the parameters of the model using the EM algorithm.
Timer::Start("em");
src/mlpack/methods/rann/ra_util.cpp
+1 -8
View File
@@ -30,11 +30,6 @@ size_t mlpack::neighbor::RAUtil::MinimumSamplesReqd(const size_t n,
double prob;
Log::Assert(alpha <= 1.0);
// going through all values of sample sizes
// to find the minimum samples required to satisfy the
// desired bound
bool done = false;
// This performs a binary search on the integer values between 'lb = k'
// and 'ub = n' to find the minimum number of samples 'm' required to obtain
// the desired success probability 'alpha'.
@@ -46,7 +41,6 @@ size_t mlpack::neighbor::RAUtil::MinimumSamplesReqd(const size_t n,
{
if (prob - alpha < 0.001 || ub < lb + 2)
{
done = true;
break;
}
else
@@ -66,12 +60,11 @@ size_t mlpack::neighbor::RAUtil::MinimumSamplesReqd(const size_t n,
}
else
{
done = true;
break;
}
}
m = (ub + lb) / 2;
} while (!done);
} while (true);
return (std::min(m + 1, n));
}
src/mlpack/tests/gmm_test.cpp
+100
View File
@@ -756,5 +756,105 @@ BOOST_AUTO_TEST_CASE(UseExistingModelTest)
}
}
/**
* Make sure we can fit a diagonal GMM reasonably.
*/
BOOST_AUTO_TEST_CASE(DiagonalGMMTrainTest)
{
Log::Warn.ignoreInput = false;
// We'll have three diagonal-covariance Gaussian distributions from this
// mixture.
distribution::GaussianDistribution d1("0.0 1.0 0.0", "1.0 0.0 0.0;"
"0.0 0.8 0.0;"
"0.0 0.0 1.0");
distribution::GaussianDistribution d2("2.0 -1.0 5.0", "3.0 0.0 0.0;"
"0.0 1.2 0.0;"
"0.0 0.0 1.3");
distribution::GaussianDistribution d3("0.0 5.0 -3.0", "2.0 0.0 0.0;"
"0.0 0.3 0.0;"
"0.0 0.0 1.0");
// Now we'll generate points and probabilities. 1500 points. Slower than I
// would like...
arma::mat points(3, 5000);
for (size_t i = 0; i < 5000; i++)
{
double randValue = math::Random();
if (randValue <= 0.20) // p(d1) = 0.20
points.col(i) = d1.Random();
else if (randValue <= 0.50) // p(d2) = 0.30
points.col(i) = d2.Random();
else // p(d3) = 0.50
points.col(i) = d3.Random();
}
// Now train the model. 3 dimensions, 3 components.
GMM g(3, 3);
g.Train<EMFit<kmeans::KMeans<>, DiagonalConstraint>>(points, 5);
// Now check the results. We need to order by weights so that when we do the
// checking, things will be correct.
arma::uvec sortedIndices = sort_index(g.Weights());
// First Gaussian (d1).
BOOST_REQUIRE_SMALL(g.Weights()[sortedIndices[0]] - 0.2, 0.1);
for (size_t i = 0; i < 3; i++)
BOOST_REQUIRE_SMALL((g.Component(sortedIndices[0]).Mean()[i]
- d1.Mean()[i]), 0.4);
for (size_t row = 0; row < 3; ++row)
{
for (size_t col = 0; col < 3; ++col)
{
const double v = g.Component(sortedIndices[0]).Covariance()(row, col);
if (row == col)
BOOST_REQUIRE_SMALL(v - d1.Covariance()(row, col), 0.5);
else
BOOST_REQUIRE_SMALL(v, 1e-5);
}
}
// Second Gaussian (d2).
BOOST_REQUIRE_SMALL(g.Weights()[sortedIndices[1]] - 0.3, 0.1);
for (size_t i = 0; i < 3; i++)
BOOST_REQUIRE_SMALL((g.Component(sortedIndices[1]).Mean()[i]
- d2.Mean()[i]), 0.4);
for (size_t row = 0; row < 3; ++row)
{
for (size_t col = 0; col < 3; ++col)
{
const double v = g.Component(sortedIndices[1]).Covariance()(row, col);
if (row == col)
BOOST_REQUIRE_SMALL(v - d2.Covariance()(row, col), 0.5);
else
BOOST_REQUIRE_SMALL(v, 1e-5);
}
}
// Third Gaussian (d3).
BOOST_REQUIRE_SMALL(g.Weights()[sortedIndices[2]] - 0.5, 0.1);
for (size_t i = 0; i < 3; ++i)
BOOST_REQUIRE_SMALL((g.Component(sortedIndices[2]).Mean()[i]
- d3.Mean()[i]), 0.4);
for (size_t row = 0; row < 3; ++row)
{
for (size_t col = 0; col < 3; ++col)
{
const double v = g.Component(sortedIndices[2]).Covariance()(row, col);
if (row == col)
BOOST_REQUIRE_SMALL(v - d3.Covariance()(row, col), 0.5);
else
BOOST_REQUIRE_SMALL(v, 1e-5);
}
}
}
BOOST_AUTO_TEST_SUITE_END();