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

Getting there.

Browse Source
This commit is contained in:
Dongryeol Lee
2010-08-16 19:38:41 +00:00
parent 3450aad895
commit 423dc84cb8
8 changed files with 182 additions and 107 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/trunk/contrib/dongryel/gp_regression/block_coord_descent.h
+38 -38
View File
@@ -69,7 +69,7 @@ class BlockCoordDescentResult: boost::noncopyable {
void Init(bool initialize_result,
const fl::data::MonolithicPoint<double> &right_hand) {
if(initialize_result) {
if (initialize_result) {
gradient_.Init(right_hand.length());
solution_.Init(right_hand.length());
}
@@ -156,7 +156,7 @@ class BlockCoordDescentDelta: boost::noncopyable {
}
template < typename Matrix, typename KernelType,
typename double >
typename double >
void Update(const Matrix &table, const KernelType &kernel,
int new_index, double self_kernel_value,
const fl::data::MonolithicPoint<double> &gradient) {
@@ -166,7 +166,7 @@ class BlockCoordDescentDelta: boost::noncopyable {
beta.Init(active_set_.size());
beta.SetZero();
for(int i = 0; i < active_set_.size(); i++) {
for (int i = 0; i < active_set_.size(); i++) {
// Compute the kernel value with the current active point.
int active_point_id = active_set_[i];
@@ -182,7 +182,7 @@ class BlockCoordDescentDelta: boost::noncopyable {
// Compute the beta by multiplying the old inverse by
// the the column we are adding into the kernel matrix.
for(int j = 0; j < active_set_.size(); j++) {
for (int j = 0; j < active_set_.size(); j++) {
beta[j] += kernel_matrix_.get(i, active_set_.size()) *
inverse_.get(j, i);
}
@@ -193,23 +193,23 @@ class BlockCoordDescentDelta: boost::noncopyable {
self_kernel_value);
// Compute eta.
for(int i = 0; i < active_set_.size(); i++) {
for (int i = 0; i < active_set_.size(); i++) {
eta += kernel_matrix_.get(active_set_.size(), i) * beta[i];
}
eta = 1.0 / (self_kernel_value - eta);
// Update the inverse.
for(int j = 0; j < active_set_.size() + 1; j++) {
for(int i = 0; i < active_set_.size() + 1; i++) {
for (int j = 0; j < active_set_.size() + 1; j++) {
for (int i = 0; i < active_set_.size() + 1; i++) {
double increment = eta;
if(i < active_set_.size()) {
if (i < active_set_.size()) {
increment *= beta[i];
}
else {
increment *= (-1);
}
if(j < active_set_.size()) {
if (j < active_set_.size()) {
increment *= beta[j];
}
else {
@@ -223,12 +223,12 @@ class BlockCoordDescentDelta: boost::noncopyable {
// Update the delta solution by computing the dot product
// between beta and the current gradient.
double factor = 0;
for(int i = 0; i < active_set_.size(); i++) {
for (int i = 0; i < active_set_.size(); i++) {
factor += beta[i] * gradient[active_set_[i]];
}
factor -= gradient[new_index];
factor *= eta;
for(int i = 0; i < active_set_.size(); i++) {
for (int i = 0; i < active_set_.size(); i++) {
delta_solution_[active_set_[i]] -= factor * beta[i];
}
delta_solution_[new_index] += factor;
@@ -252,7 +252,7 @@ class BlockCoordDescentDelta: boost::noncopyable {
temp_gradient_info_.SetZero();
inactive_set_.resize(num_points_);
active_set_.resize(0);
for(int i = 0; i < num_points_; i++) {
for (int i = 0; i < num_points_; i++) {
inactive_set_[i] = i;
}
inverse_.SetZero();
@@ -292,9 +292,9 @@ class SolveSubProblem {
const std::vector<int> &active_set = delta.active_set();
fl::dense::Matrix<double, false> &kernel_matrix = delta.kernel_matrix();
for(int j = 0; j < active_set.size(); j++) {
for (int j = 0; j < active_set.size(); j++) {
int column_index = active_set[j];
for(int i = 0; i < active_set.size(); i++) {
for (int i = 0; i < active_set.size(); i++) {
int row_index = active_set[i];
kernel_matrix.set(
i, j, fl::ml::KernelValue::Compute(
@@ -310,10 +310,10 @@ class SolveSubProblem {
fl::data::MonolithicPoint<double> residual;
residual.Alias(delta.temp_gradient_info().ptr(), kernel_matrix.n_rows());
for(int i = 0; i < active_set.size(); i++) {
for (int i = 0; i < active_set.size(); i++) {
int row_index = active_set[i];
residual[i] = right_hand_side[row_index];
for(int j = 0; j < active_set.size(); j++) {
for (int j = 0; j < active_set.size(); j++) {
int col_index = active_set[j];
residual[i] -= kernel_matrix.get(i, j) * current_solution[col_index];
}
@@ -326,7 +326,7 @@ class SolveSubProblem {
fl::data::MonolithicPoint<double> &delta_solution =
delta.delta_solution();
delta_solution.SetZero();
for(int i = 0; i < kernel_matrix.n_rows(); i++) {
for (int i = 0; i < kernel_matrix.n_rows(); i++) {
int row_index = active_set[i];
delta_solution[row_index] = -tmp_vector_[i];
}
@@ -365,8 +365,8 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GREEDY_BC> {
max_random_set_size,
(int)inactive_set.size());
for(int i = ((int) inactive_set.size()) - 1;
i >= ((int) inactive_set.size()) - (*random_set_size); i--) {
for (int i = ((int) inactive_set.size()) - 1;
i >= ((int) inactive_set.size()) - (*random_set_size); i--) {
int random_index = fl::math::Random(0, i + 1);
std::swap(inactive_set[random_index], inactive_set[i]);
@@ -426,22 +426,22 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GREEDY_BC> {
double minimum_value =
std::numeric_limits<double>::max();
double cached_self_kernel_value = -1;
for(int i = inactive_set.size() - random_set_size;
i < inactive_set.size(); i++) {
for (int i = inactive_set.size() - random_set_size;
i < inactive_set.size(); i++) {
double kernel_value =
fl::ml::KernelValue::Compute(
table, kernel, inactive_set[i], inactive_set[i], true);
double current_value =
- fl::math::Sqr(temp_gradient_info[inactive_set[i]]) /
(2.0 * kernel_value);
if(current_value <= minimum_value) {
if (current_value <= minimum_value) {
minimum_value = current_value;
selected_index = i;
cached_self_kernel_value = kernel_value;
}
}
if(delta.current_active_set_size() == 0) {
if (delta.current_active_set_size() == 0) {
inverse.set(0, 0, 1.0 / cached_self_kernel_value);
kernel_matrix.set(0, 0, cached_self_kernel_value);
delta_solution[inactive_set[selected_index]] =
@@ -466,15 +466,15 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GREEDY_BC> {
ChooseSubset(
result, inactive_set, max_random_set_size, &random_set_size);
for(int i = inactive_set.size() - random_set_size;
i < inactive_set.size(); i++) {
for (int i = inactive_set.size() - random_set_size;
i < inactive_set.size(); i++) {
// The index of the randomly chosen point.
int random_point_index = inactive_set[i];
double dot_product = 0;
// Loop over each point in the current active set.
for(int j = 0; j < active_set.size(); j++) {
for (int j = 0; j < active_set.size(); j++) {
int active_point_index = active_set[j];
double kernel_value =
fl::ml::KernelValue::Compute(
@@ -487,8 +487,8 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GREEDY_BC> {
}
}
while(delta.current_active_set_size() < active_set_size &&
inactive_set.size() > 0);
while (delta.current_active_set_size() < active_set_size &&
inactive_set.size() > 0);
}
};
@@ -521,12 +521,12 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::CYCLIC_BC> {
std::vector<int> &inactive_set = delta.inactive_set();
active_set.resize(0);
inactive_set.resize(0);
for(int i = 0; i < std::min(active_set_size, table.n_entries()); i++) {
for (int i = 0; i < std::min(active_set_size, table.n_entries()); i++) {
int active_index = (i + starting_index_) % table.n_entries();
active_set.push_back(active_index);
}
for(int i = std::min(active_set_size, table.n_entries());
i < table.n_entries(); i++) {
for (int i = std::min(active_set_size, table.n_entries());
i < table.n_entries(); i++) {
int inactive_index = (i + starting_index_) % table.n_entries();
inactive_set.push_back(inactive_index);
}
@@ -591,7 +591,7 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GRADIENT_BC> {
// Sort the current gradient component by absolute magnitude and
// select the indices with the greatest ones.
for(int i = 0; i < sorted_indices_.size(); i++) {
for (int i = 0; i < sorted_indices_.size(); i++) {
sorted_indices_[i] = i;
}
@@ -603,11 +603,11 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GRADIENT_BC> {
std::vector<int> &inactive_set = delta.inactive_set();
active_set.resize(0);
inactive_set.resize(0);
for(int i = 0; i < std::min(active_set_size, table.n_entries()); i++) {
for (int i = 0; i < std::min(active_set_size, table.n_entries()); i++) {
active_set.push_back(sorted_indices_[i]);
}
for(int i = std::min(active_set_size, table.n_entries());
i < table.n_entries(); i++) {
for (int i = std::min(active_set_size, table.n_entries());
i < table.n_entries(); i++) {
inactive_set.push_back(sorted_indices_[i]);
}
@@ -639,11 +639,11 @@ class BlockCoordDescent: private boost::noncopyable {
fl::data::MonolithicPoint<double> &solution = result.solution();
fl::data::MonolithicPoint<double> &gradient = result.gradient();
for(int i = 0; i < active_set.size(); i++) {
for (int i = 0; i < active_set.size(); i++) {
solution[active_set[i]] -= delta_solution[active_set[i]];
double dot_product = 0;
for(int j = 0; j < active_set.size(); j++) {
for (int j = 0; j < active_set.size(); j++) {
dot_product += kernel_matrix.get(i, j) *
delta_solution[active_set[j]];
}
@@ -674,7 +674,7 @@ class BlockCoordDescent: private boost::noncopyable {
active_set_size, table.n_entries());
// The main loop.
while(!result->IsConverged()) {
while (!result->IsConverged()) {
active_set_selection_trait.Select(
active_set_size, max_random_set_size, table, kernel, right_hand,
fastlib/trunk/contrib/dongryel/gp_regression/dense_matrix_inverse.h
+3 -3
View File
@@ -24,8 +24,8 @@ class DenseMatrixInverse {
new_matrix_inverse->Init(previous_inverse.n_rows() + 1,
previous_inverse.n_cols() + 1);
for(int j = 0; j < previous_inverse.n_cols(); j++) {
for(int i = 0; i < previous_inverse.n_rows(); i++) {
for (int j = 0; j < previous_inverse.n_cols(); j++) {
for (int i = 0; i < previous_inverse.n_rows(); i++) {
new_matrix_inverse->set(
i, j, previous_inverse.get(i, j) +
inverse_times_new_column[i] *
@@ -33,7 +33,7 @@ class DenseMatrixInverse {
}
}
for(int j = 0; j < previous_inverse.n_cols(); j++) {
for (int j = 0; j < previous_inverse.n_cols(); j++) {
new_matrix_inverse->set(
j,
previous_inverse.n_cols(), - inverse_times_new_column[j] /
fastlib/trunk/contrib/dongryel/gp_regression/dictionary_dev.h
+11 -11
View File
@@ -28,10 +28,10 @@ Dictionary::Dictionary() {
}
Dictionary::~Dictionary() {
if(current_kernel_matrix_ != NULL) {
if (current_kernel_matrix_ != NULL) {
delete current_kernel_matrix_;
}
if(current_kernel_matrix_inverse_ != NULL) {
if (current_kernel_matrix_inverse_ != NULL) {
delete current_kernel_matrix_inverse_;
}
}
@@ -57,8 +57,8 @@ void Dictionary::inactive_indices(
// Scan the in_dictionary list and build the inactive index set.
inactive_indices_out->resize(0);
for(int i = 0; i < in_dictionary_.size(); i++) {
if(in_dictionary_[i] == false) {
for (int i = 0; i < in_dictionary_.size(); i++) {
if (in_dictionary_[i] == false) {
inactive_indices_out->push_back(i);
}
}
@@ -83,12 +83,12 @@ void Dictionary::UpdateDictionary_(
current_kernel_matrix_->n_rows() + 1,
current_kernel_matrix_->n_cols() + 1);
for(int j = 0; j < current_kernel_matrix_->n_cols(); j++) {
for(int i = 0; i < current_kernel_matrix_->n_rows(); i++) {
for (int j = 0; j < current_kernel_matrix_->n_cols(); j++) {
for (int i = 0; i < current_kernel_matrix_->n_rows(); i++) {
new_kernel_matrix->set(i, j, current_kernel_matrix_->get(i, j));
}
}
for(int j = 0; j < current_kernel_matrix_->n_cols(); j++) {
for (int j = 0; j < current_kernel_matrix_->n_cols(); j++) {
new_kernel_matrix->set(
j, current_kernel_matrix_->n_cols(), new_column_vector[j]);
new_kernel_matrix->set(
@@ -138,11 +138,11 @@ void Dictionary::AddBasis(
const std::vector<double> &new_column_vector_in,
double self_value) {
if(new_column_vector_in.size() > 0) {
if (new_column_vector_in.size() > 0) {
Vector new_column_vector;
new_column_vector.Init(new_column_vector_in.size());
for(int i = 0; i < new_column_vector.length(); i++) {
for (int i = 0; i < new_column_vector.length(); i++) {
new_column_vector[i] = new_column_vector_in[i];
}
@@ -160,7 +160,7 @@ void Dictionary::AddBasis(
// If the projection error is above the threshold, add it to the
// dictionary.
if(projection_error > adding_threshold_) {
if (projection_error > adding_threshold_) {
UpdateDictionary_(
new_point_index,
new_column_vector,
@@ -191,7 +191,7 @@ void Dictionary::Init(const Matrix *table_in) {
// Generate a random permutation and initialize the inital
// dictionary which consists of the first random point.
random_permutation_.resize(table_in->n_cols());
for(int i = 0; i < table_in->n_cols(); i++) {
for (int i = 0; i < table_in->n_cols(); i++) {
random_permutation_[i] = i;
in_dictionary_[i] = false;
training_index_to_dictionary_position_[i] = -1;
fastlib/trunk/contrib/dongryel/gp_regression/kernel_linear_operator.h
+10 -10
View File
@@ -45,7 +45,7 @@ class DotProductTrait<false> {
};
template < typename TableType, typename KernelType, bool do_centering,
bool dotproduct_selfcase_special = false >
bool dotproduct_selfcase_special = false >
class KernelLinearOperator: public LinearOperator {
private:
@@ -74,23 +74,23 @@ class KernelLinearOperator: public LinearOperator {
comm_ = &comm_in;
map_ = &map_in;
if(do_centering) {
if (do_centering) {
average_row_.Init(table_in.n_entries());
}
average_ = 0;
if(do_centering) {
if (do_centering) {
// Precompute the average. This is a naive way of computing it.
for(int i = 0; i < table_in.n_entries(); i++) {
for (int i = 0; i < table_in.n_entries(); i++) {
double average_for_i_th_point = 0;
for(int j = 0; j < table_in.n_entries(); j++) {
for (int j = 0; j < table_in.n_entries(); j++) {
average_for_i_th_point += kernel_value(i, j);
}
average_for_i_th_point /= ((double) table_in.n_entries());
average_row_[i] = average_for_i_th_point;
}
for(int i = 0; i < table_in.n_entries(); i++) {
for (int i = 0; i < table_in.n_entries(); i++) {
average_ += average_row_[i];
}
average_ /= ((double) table_in.n_entries());
@@ -125,11 +125,11 @@ class KernelLinearOperator: public LinearOperator {
prods.PutScalar(0);
for(int j = 0; j < table_->n_entries(); j++) {
for(int i = 0; i < table_->n_entries(); i++) {
for (int j = 0; j < table_->n_entries(); j++) {
for (int i = 0; i < table_->n_entries(); i++) {
double pair_kernel_value = (do_centering) ?
centered_kernel_value(i, j) : kernel_value(i, j);
for(int k = 0; k < vecs.NumVectors(); k++) {
for (int k = 0; k < vecs.NumVectors(); k++) {
prods.Pointers()[k][i] += pair_kernel_value * vecs.Pointers()[k][j];
}
}
@@ -152,7 +152,7 @@ class KernelLinearOperator: public LinearOperator {
template<typename TableType, typename KernelType, bool do_centering>
class OperatorTraits < double, Epetra_MultiVector,
KernelLinearOperator<TableType, KernelType, do_centering> > {
KernelLinearOperator<TableType, KernelType, do_centering> > {
public:
static void Apply(const Epetra_Operator& Op,
fastlib/trunk/contrib/dongryel/gp_regression/lbfgs_dev.h
+16 -16
View File
@@ -60,7 +60,7 @@ double Lbfgs<FunctionType>::ChooseScalingFactor_(
const Vector &gradient) {
double scaling_factor = 1.0;
if(iteration_num > 0) {
if (iteration_num > 0) {
int previous_pos = (iteration_num - 1) % num_basis_;
Vector s_basis;
Vector y_basis;
@@ -113,7 +113,7 @@ bool Lbfgs<FunctionType>::LineSearch_(
la::Dot(gradient, search_direction);
// If it is not a descent direction, just report failure.
if(initial_search_direction_dot_gradient > 0.0) {
if (initial_search_direction_dot_gradient > 0.0) {
return false;
}
@@ -131,7 +131,7 @@ bool Lbfgs<FunctionType>::LineSearch_(
const double inc = 2.1;
const double dec = 0.5;
double width = 0;
for(; ;) {
for (; ;) {
// Perform a step and evaluate the gradient and the function
// values at that point.
@@ -141,7 +141,7 @@ bool Lbfgs<FunctionType>::LineSearch_(
function_->Gradient(new_iterate_tmp_, &gradient);
num_iterations++;
if(function_value > initial_function_value + step_size *
if (function_value > initial_function_value + step_size *
linear_approx_function_value_decrease) {
width = dec;
}
@@ -151,12 +151,12 @@ bool Lbfgs<FunctionType>::LineSearch_(
double search_direction_dot_gradient =
la::Dot(gradient, search_direction);
if(search_direction_dot_gradient < param_.wolfe() *
if (search_direction_dot_gradient < param_.wolfe() *
initial_search_direction_dot_gradient) {
width = inc;
}
else {
if(search_direction_dot_gradient > -param_.wolfe() *
if (search_direction_dot_gradient > -param_.wolfe() *
initial_search_direction_dot_gradient) {
width = dec;
}
@@ -168,13 +168,13 @@ bool Lbfgs<FunctionType>::LineSearch_(
// Terminate when the step size gets too small or too big or it
// exceeds the max number of iterations.
if(step_size < param_.min_step()) {
if (step_size < param_.min_step()) {
return false;
}
if(step_size > param_.max_step()) {
if (step_size > param_.max_step()) {
return false;
}
if(num_iterations >= param_.max_line_search()) {
if (num_iterations >= param_.max_line_search()) {
return false;
}
@@ -204,7 +204,7 @@ void Lbfgs<FunctionType>::SearchDirection_(
int limit = std::max(iteration_num - num_basis_, 0);
for(int i = iteration_num - 1; i >= limit; i--) {
for (int i = iteration_num - 1; i >= limit; i--) {
int translated_position = i % num_basis_;
Vector y_basis, s_basis;
s_lbfgs_.MakeColumnVector(translated_position, &s_basis);
@@ -214,7 +214,7 @@ void Lbfgs<FunctionType>::SearchDirection_(
la::Dot(s_basis, q);
}
la::ScaleOverwrite(scaling_factor, q, search_direction);
for(int i = limit; i <= iteration_num - 1; i++) {
for (int i = limit; i <= iteration_num - 1; i++) {
int translated_position = i % num_basis_;
Vector y_basis, s_basis;
s_lbfgs_.MakeColumnVector(translated_position, &s_basis);
@@ -254,7 +254,7 @@ double Lbfgs<FunctionType>::Evaluate_(
// value encountered during the optimization.
double function_value = function_->Evaluate(iterate);
if(function_value < min_point_iterate_.second) {
if (function_value < min_point_iterate_.second) {
min_point_iterate_.first.CopyValues(iterate);
min_point_iterate_.second = function_value;
}
@@ -311,11 +311,11 @@ bool Lbfgs<FunctionType>::Optimize(int num_iterations,
// The main optimization loop.
int it_num;
for(it_num = 0; optimize_until_convergence ||
it_num < num_iterations; it_num++) {
for (it_num = 0; optimize_until_convergence ||
it_num < num_iterations; it_num++) {
// Break when the norm of the gradient becomes too small.
if(GradientNormTooSmall_(gradient)) {
if (GradientNormTooSmall_(gradient)) {
break;
}
@@ -336,7 +336,7 @@ bool Lbfgs<FunctionType>::Optimize(int num_iterations,
LineSearch_(function_value, *iterate, gradient, search_direction,
step_size);
if(search_is_success == false) {
if (search_is_success == false) {
break;
}
fastlib/trunk/contrib/dongryel/gp_regression/linear_operator.h
+2 -2
View File
@@ -110,8 +110,8 @@ class LinearOperator: public virtual Epetra_Operator {
void PrintDebug(const char *name = "", FILE *stream = stderr) const {
fprintf(stream, "----- MATRIX ------: %s\n", name);
for(int r = 0; r < this->n_rows(); r++) {
for(int c = 0; c < this->n_cols(); c++) {
for (int r = 0; r < this->n_rows(); r++) {
for (int c = 0; c < this->n_cols(); c++) {
fprintf(stream, "%+3.3f ", this->get(r, c));
}
fprintf(stream, "\n");
fastlib/trunk/contrib/dongryel/gp_regression/sg_gp_regression.h
+15 -2
View File
@@ -38,14 +38,27 @@ class SparseGreedyGprModel {
const Vector &right_hand_side_in,
Vector *solution_out) const;
void ExtractWeightedTargetSubset_(
const ml::Dictioanay &dictionary_in,
const std::vector<double> &additional_kernel_values_in,
Vector *target_subset_out) const;
void ExtractTargetSubset_(
const ml::Dictioanay &dictionary_in,
Vector *target_subset_out) const;
double QuadraticObjective_(
const ml::Dictionary &dictionary_in,
const std::vector<double> &kernel_values_in,
bool for_coeffs) const;
template<typename CovarianceType>
void ComputeKernelValues_(
const CovarianceType &covariance_in,
const std::vector<int> &point_indices_in_dictionary,
const Vector &point,
Vector *kernel_values_out) const;
template<typename CovarianceType>
void ComputeKernelValues_(
const CovarianceType &covariance_in,
@@ -84,12 +97,12 @@ class SparseGreedyGprModel {
void FinalizeModel();
template<typename CovarianceType>
void PredictMean(
double PredictMean(
const CovarianceType &covariance,
const Vector &point) const;
template<typename CovarianceType>
void PredictVariance(
double PredictVariance(
const CovarianceType &covariance,
const Vector &point) const;
};
fastlib/trunk/contrib/dongryel/gp_regression/sg_gp_regression_dev.h
+87 -25
View File
@@ -15,27 +15,45 @@ namespace ml {
namespace gp_regression {
template<typename CovarianceType>
void SparseGreedyGprModel::PredictMean(
double SparseGreedyGprModel::PredictMean(
const CovarianceType &covariance,
const Vector &point) const {
// Get the indices of the points in the dictionary.
const std::vector<int> &point_indices_in_dictionary =
dictionary_for_error_.point_indices_in_dictionary();
// Compute the kernel values with the basis points and take the dot
// product with the coefficients.
Vector kernel_values;
ComputeKernelValues_(
covariance, point_indices_in_dictionary, point, &kernel_values);
return la::Dot(kernel_values, coefficients_);
}
template<typename CovarianceType>
void SparseGreedyGprModel::PredictVariance(
double SparseGreedyGprModel::PredictVariance(
const CovarianceType &covariance,
const Vector &point) const {
// Get the indices of the points in the dictionary.
const std::vector<int> &point_indices_in_dictionary =
dictionary_for_error_.point_indices_in_dictionary();
// The variance is self-correlation minus the variance explained by
// the sparse GPR model.
Vector kernel_values;
ComputeKernelValues_(
covariance, point_indices_in_dictionary, point, &kernel_values);
}
void SparseGreedyGprModel::FinalizeModel() {
SolveSystem_(dictionary_, &coefficients_);
Vector weighted_target_subset;
ExtractWeightedTargetSubset_(
dictionary_in, kernel_values_in, &weighted_target_subset);
SolveSystem_(dictionary_, weighted_target_subset, &coefficients);
}
double SparseGreedyGprModel::frobenius_norm_targets() const {
@@ -53,6 +71,30 @@ void SparseGreedyGprModel::SolveSystem_(
solution_out);
}
void SparseGreedyGprModel::ExtractWeightedTargetSubset_(
const ml::Dictioanay &dictionary_in,
const std::vector<double> &additional_kernel_values_in,
Vector *target_subset_out) const {
// Get the indices of the points in the dictionary.
const std::vector<int> &point_indices_in_dictionary =
dictionary_in.point_indices_in_dictionary();
target_subset_out->Init(point_indices_in_dictionary.size());
for (int i = 0; i < target_subset_out->length(); i++) {
int basis_index = point_indices_in_dictionary[i];
const std::vector<double> &kernel_values =
(i == target_subset_out->length() - 1) ?
additional_kernel_values_in : kernel_matrix_columns_[basis_index];
double dot_product = 0;
for (int j = 0; j < dataset_->n_cols(); j++) {
dot_product += (*targets_)[ j ] * kernel_values[j];
}
(*target_subset_out)[i] = dot_product;
}
}
void SparseGreedyGprModel::ExtractTargetSubset_(
const ml::Dictioanay &dictionary_in,
Vector *target_subset_out) const {
@@ -62,28 +104,31 @@ void SparseGreedyGprModel::ExtractTargetSubset_(
dictionary_in.point_indices_in_dictionary();
target_subset_out->Init(point_indices_in_dictionary.size());
for(int i = 0; i < target_subset_out->length(); i++) {
for (int i = 0; i < target_subset_out->length(); i++) {
(*target_subset_out)[i] = (*targets_)[ point_indices_in_dictionary[i] ];
}
}
double SparseGreedyGprModel::QuadraticObjective_(
const ml::Dictionary &dictionary_in,
const std::vector<double> &kernel_values_in,
bool for_coeffs) const {
// Compute the objective, -0.5 y^T K^{1} y.
Vector product;
if(for_coeffs) {
if (for_coeffs) {
Vector weighted_target_subset;
ExtractWeightedTargetSubset_(
dictionary_in, kernel_values_in, &weighted_target_subset);
SolveSystem_(dictionary_in, weighted_target_subset, &product);
return -0.5 * la::Dot(product, weighted_target_subset);
}
else {
Vector target_subset;
ExtractTargetSubset_(dictionary_in, &target_subset);
SolveSystem_(dictionary_in, target_subset, &product);
return -0.5 * la::Dot(product, target_subset);
}
return - 0.5 * la::Dot(product, target_subset);
}
void SparseGreedyGprModel::FillSquaredKernelMatrix_(
@@ -97,15 +142,15 @@ void SparseGreedyGprModel::FillSquaredKernelMatrix_(
const std::vector<int> &point_indices_in_dictionary =
dictionary_.point_indices_in_dictionary();
for(int i = 0; i < point_indices_in_dictionary.size(); i++) {
for (int i = 0; i < point_indices_in_dictionary.size(); i++) {
int basis_index = point_indices_in_dictionary[i];
const std::vector<double> &cached_kernel_values =
kernel_matrix_columns_[basis_index];
// Take the dot product.
double dot_product = 0;
for(int j = 0; j < kernel_values.length(); j++) {
if(j == candidate_index) {
for (int j = 0; j < kernel_values.length(); j++) {
if (j == candidate_index) {
dot_product += cached_kernel_values[j] *
(kernel_values[j] + noise_level_in);
}
@@ -132,13 +177,29 @@ void SparseGreedyGprModel::FillKernelMatrix_(
// Simply the necessary kernel values. For the self value, add the
// noise.
for(int i = 0; i < point_indices_in_dictionary.size(); i++) {
for (int i = 0; i < point_indices_in_dictionary.size(); i++) {
(*new_column_vector_out)[i] =
kernel_values[ point_indices_in_dictionary[i] ];
}
*new_self_value_out = kernel_values[candidate_index] + noise_level_in;
}
template<typename CovarianceType>
void SparseGreedyGprModel::ComputeKernelValues_(
const CovarianceType &covariance_in,
const std::vector<int> &point_indices_in_dictionary,
const Vector &point,
Vector *kernel_values_out) const {
kernel_values_out->Init(point_indices_in_dictionary.size());
for (int i = 0; i < kernel_values_out->length(); i++) {
int basis_index = point_indices_in_dictionary[i];
Vector basis_point;
dataset_->MakeColumnVector(i, &basis_point);
(*kernel_values_out)[i] = covariance_in.Dot(point, basis_point, false);
}
}
template<typename CovarianceType>
void SparseGreedyGprModel::ComputeKernelValues_(
const CovarianceType &covariance_in,
@@ -149,7 +210,7 @@ void SparseGreedyGprModel::ComputeKernelValues_(
dataset_->MakeColumnVector(candidate_index, &candidate_point);
// Fill out the kernel values sequentially.
for(int i = 0; i < dataset_->n_cols(); i++) {
for (int i = 0; i < dataset_->n_cols(); i++) {
Vector point;
dataset_->MakeColumnVector(i, &point);
(*kernel_values_out)[i] = covariance_in.Dot(
@@ -172,11 +233,11 @@ double SparseGreedyGprModel::AddOptimalPoint(
double optimum_value = std::numeric_limits<double>::max();
// Loop over candidates and decide to add the optimal.
for(int i = 0; i < candidate_indices.size(); i++) {
for (int i = 0; i < candidate_indices.size(); i++) {
// Make a copy of the dictionaries.
Dictionary dictionary_copy;
if(for_coeffs) {
if (for_coeffs) {
dictionary_copy = dictionary_;
}
else {
@@ -193,7 +254,7 @@ double SparseGreedyGprModel::AddOptimalPoint(
// Compute the additional quantities to be appended to grow the
// matrix and update the dictionary.
if(for_coeffs) {
if (for_coeffs) {
FillSquaredKernelMatrix_(
candidate_index, kernel_values, &new_column_vector, &new_self_value);
}
@@ -205,9 +266,10 @@ double SparseGreedyGprModel::AddOptimalPoint(
candidate_index, new_column_vector, new_self_value);
// Compute the objective function value for the coefficients.
double objective_value = QuadraticObjective_(dictionary_copy, for_coeffs);
double objective_value = QuadraticObjective_(
dictionary_copy, kernel_values, for_coeffs);
if(objective_value < optimum_value) {
if (objective_value < optimum_value) {
optimal_point_index = i;
optimum_value = objective_value;
}
@@ -217,7 +279,7 @@ double SparseGreedyGprModel::AddOptimalPoint(
int final_candidate_index = candidate_indices[optimal_point_index];
std::vector<double> final_column_vector;
double final_self_value;
if(for_coeffs) {
if (for_coeffs) {
// Grow the kernel cache in this case.
kernel_matrix_columns_[final_candidate_index].resize(dataset_->n_cols());
@@ -265,13 +327,13 @@ void SparseGreedyGpr::ChooseRandomSubset_(
// Copy
subset_out->resize(inactive_set.size());
for(int i = 0; i < inactive_set.size(); i++) {
for (int i = 0; i < inactive_set.size(); i++) {
(*subset_out)[i] = inactive_set[i];
}
// Then shuffle, and truncate.
if(inactive_set.size() > subset_size) {
for(int i = subset_out->size() - 1; i >= 1; i--) {
if (inactive_set.size() > subset_size) {
for (int i = subset_out->size() - 1; i >= 1; i--) {
// Pick a random index between 0 and i, inclusive.
int random_index = math::RandInt(0, i);
@@ -321,7 +383,7 @@ void SparseGreedyGpr::Compute(
// sets).
std::vector<int> inactive_indices(dataset_->n_cols());
std::vector<int> inactive_indices_for_error(dataset_-> n_cols());
for(int i = 0; i < dataset_->n_cols(); i++) {
for (int i = 0; i < dataset_->n_cols(); i++) {
inactive_indices[i] = i;
inactive_indices_for_error[i] = i;
}
@@ -351,9 +413,9 @@ void SparseGreedyGpr::Compute(
dictionary_.inactive_indices(&inactive_indices);
dictionary_for_error_.inactive_indices(&inactive_indices_for_error_);
}
while(Done_(
model_out->frobenius_norm_targets(), noise_level_in, precision_in,
optimum_value, optimum_value_for_error));
while (Done_(
model_out->frobenius_norm_targets(), noise_level_in, precision_in,
optimum_value, optimum_value_for_error));
// Using the final model, compute the coefficients.
model_out->FinalizeModel();