From fb5adf58b79f7ec4bc16ec489cbba8f9df8e5e7f Mon Sep 17 00:00:00 2001 From: Dongryeol Lee Date: Sun, 24 Feb 2008 21:52:16 +0000 Subject: [PATCH] Removed the custom-built Krylov solver, now going to replace it with Epetra solver --- fastlib2/contrib/dongryel/regression/build.py | 1 - .../contrib/dongryel/regression/krylov_lpr.h | 225 ++--- .../regression/krylov_lpr_setup_impl.h | 86 +- .../regression/krylov_lpr_solver_impl.h | 857 ------------------ .../dongryel/regression/krylov_lpr_test.h | 72 +- .../contrib/dongryel/regression/krylov_stat.h | 100 -- .../dongryel/regression/relative_prune_lpr.h | 88 +- 7 files changed, 121 insertions(+), 1308 deletions(-) delete mode 100644 fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h diff --git a/fastlib2/contrib/dongryel/regression/build.py b/fastlib2/contrib/dongryel/regression/build.py index 70686aba3f..250da2bc1f 100644 --- a/fastlib2/contrib/dongryel/regression/build.py +++ b/fastlib2/contrib/dongryel/regression/build.py @@ -21,7 +21,6 @@ librule( headers = ["epan_kernel_moment_info.h", "krylov_lpr.h", "krylov_lpr_setup_impl.h", - "krylov_lpr_solver_impl.h", "krylov_lpr_test.h", "lpr_util.h", "naive_lpr.h"], diff --git a/fastlib2/contrib/dongryel/regression/krylov_lpr.h b/fastlib2/contrib/dongryel/regression/krylov_lpr.h index 2be89f8dde..582d1df878 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr.h +++ b/fastlib2/contrib/dongryel/regression/krylov_lpr.h @@ -16,6 +16,7 @@ #include "epan_kernel_moment_info.h" #include "multi_index_util.h" #include "lpr_util.h" +#include "mlpack/allknn/allknn.h" #define INSIDE_KRYLOV_LPR_H #include "krylov_stat.h" @@ -72,12 +73,23 @@ class KrylovLpr { */ Vector rset_targets_; + /** @brief The reference training target value divided by the + * normalization constant of the kernel centered at each + * reference point. + */ + Vector rset_target_divided_by_norm_consts_; + + /** @brief The inverse of the normalization constant of the kernel + * centered at each reference point. + */ + Vector rset_inv_norm_consts_; + /** @brief The original training target value for the reference * dataset weighted by the reference coordinate. (i.e. y_i * [1; r^T]^T ). */ Matrix target_weighted_rset_; - + /** @brief The computed fit values at each reference point. */ Vector rset_regression_estimates_; @@ -148,43 +160,16 @@ class KrylovLpr { void TestRightHandSideComputation_(const Matrix &qset, const Matrix &approximated); - /** @brief This function test the second phase computation (i.e. - * the computation of the product of B^T W(q) B and z(q). - */ - void TestKrylovComputation_ - (const Matrix &qset, const Matrix &approximated, - const Matrix ¤t_lanczos_vectors, - const ArrayList &query_should_exit_the_loop); - - void NormalizeMatrixColumnVectors_(Matrix &m, Vector &lengths) { - - for(index_t i = 0; i < m.n_cols(); i++) { - double *column_vector = m.GetColumnPtr(i); - lengths[i] = la::LengthEuclidean(row_length_, column_vector); - - if(lengths[i] > 0) { - la::Scale(row_length_, 1.0 / lengths[i], column_vector); - } - } - } - - /** @brief Compute the dot-product bounds possible for a pair of - * point lying in each of the two given regions. - */ - void DotProductBetweenTwoBounds_(QueryTree *qnode, ReferenceTree *rnode, - DRange &negative_dot_product_range, - DRange &positive_dot_product_range); - /** @brief Initialize the bound statistics relevant to the right * hand side computation. */ - void InitializeQueryTreeRightHandSides_(QueryTree *qnode); + void InitializeQueryTree_(QueryTree *qnode); /** @brief The postprocessing function to finalize the computation * of the right-hand sides of the linear system for each * query point. */ - void FinalizeQueryTreeRightHandSides_ + void FinalizeQueryTree_ (QueryTree *qnode, const Matrix &qset, Matrix &right_hand_sides_l, Matrix &right_hand_sides_e, Vector &right_hand_sides_used_error, Vector &right_hand_sides_n_pruned); @@ -192,16 +177,17 @@ class KrylovLpr { /** @brief Preprocess the reference tree for bottom up statistics * computation. */ - void InitializeReferenceStatistics_(ReferenceTree *rnode); + void InitializeReferenceStatistics_(ReferenceTree *rnode, int column_index, + const Vector &weights); /** @brief Determine whether the given query and the reference node * pair can be pruned. * * @return True, if it can be pruned. False, otherwise. */ - bool PrunableRightHandSides_(QueryTree *qnode, ReferenceTree *rnode, - DRange &dsqd_range, DRange &kernel_value_range, - double &used_error); + bool PrunableKrylov_(QueryTree *qnode, ReferenceTree *rnode, + DRange &dsqd_range, DRange &kernel_value_range, + double &used_error); /** @brief The base-case exhaustive computation for dual-tree based * computation of B^T W(q) Y. @@ -209,7 +195,7 @@ class KrylovLpr { * @param qnode The query node. * @param rnode The reference node. */ - void DualtreeRightHandSidesBase_ + void DualtreeWeightedVectorSumBase_ (QueryTree *qnode, ReferenceTree *rnode, const Matrix &qset, Matrix &right_hand_sides_l, Matrix &right_hand_sides_e, Vector &right_hand_sides_used_error, Vector &right_hand_sides_n_pruned); @@ -220,7 +206,7 @@ class KrylovLpr { * @param qnode The query node. * @param rnode The reference node. */ - void DualtreeRightHandSidesCanonical_ + void DualtreeWeightedVectorSumCanonical_ (QueryTree *qnode, ReferenceTree *rnode, const Matrix &qset, Matrix &right_hand_sides_l, Matrix &right_hand_sides_e, Vector &right_hand_sides_used_error, Vector &right_hand_sides_n_pruned); @@ -231,106 +217,31 @@ class KrylovLpr { * z(q) = B^T W(q) Y. This function calls a dual-tree based * fast vector summation to achieve this effect. */ - void ComputeRightHandSides_ - (QueryTree *qroot, const Matrix &qset, Matrix &right_hand_sides_l, - Matrix &right_hand_sides_e, Vector &right_hand_sides_used_error, - Vector &right_hand_sides_n_pruned) { + void ComputeWeightedVectorSum_ + (QueryTree *qroot, const Matrix &qset, const Vector &weights, + Matrix &right_hand_sides_l, Matrix &right_hand_sides_e, + Vector &right_hand_sides_used_error, Vector &right_hand_sides_n_pruned) { - // Initialize the bound quantities. + // Initialize the weight statistics on the reference side. + InitializeReferenceStatistics_(rroot_, 0, weights); + + // Initialize the bound quantities on the query side. right_hand_sides_l.SetZero(); right_hand_sides_e.SetZero(); right_hand_sides_used_error.SetZero(); right_hand_sides_n_pruned.SetZero(); - InitializeQueryTreeRightHandSides_(qroot); - + InitializeQueryTree_(qroot); + // Call dualtree function. - DualtreeRightHandSidesCanonical_ + DualtreeWeightedVectorSumCanonical_ (qroot, rroot_, qset, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); // Final traversal of the query tree to finalize estimates. - FinalizeQueryTreeRightHandSides_ - (qroot, qset, right_hand_sides_l, right_hand_sides_e, - right_hand_sides_used_error, right_hand_sides_n_pruned); + FinalizeQueryTree_(qroot, qset, right_hand_sides_l, right_hand_sides_e, + right_hand_sides_used_error, right_hand_sides_n_pruned); } - /** @brief Initialize the query tree for an iteration inside a - * Krylov solver. This forms the bounds for the solution - * vectors owned by the query points for a given query node. - * - * @param qnode The current query node. - * @param current_lanczos_vectors Each column of this matrix is a current - * Lanczos vector for each query point. - */ - void InitializeQueryTreeLanczosVectorBound_ - (QueryTree *qnode, const Matrix &qset, - const ArrayList &exclude_query_flag, - const Matrix ¤t_lanczos_vectors); - - /** @brief Finalize the Lanczos vector generator by traversing the - * query tree and summing up any unincorporated quantities. - * - * @param qnode The query node. - */ - void FinalizeQueryTreeLanczosMultiplier_ - (QueryTree *qnode, const Matrix &qset, - const ArrayList &exclude_query_flag, - const Matrix ¤t_lanczos_vectors, - Matrix &lanczos_prod_l, Matrix &lanczos_prod_e, - Vector &lanczos_prod_used_error, Vector &lanczos_prod_n_pruned, - Matrix &neg_lanczos_prod_e, Matrix &neg_lanczos_prod_u, - Vector &neg_lanczos_prod_used_error, Vector &neg_lanczos_prod_n_pruned); - - /** @brief Determine whether the given query and the reference node - * pair can be pruned. - * - * @return True, if it can be pruned. False, otherwise. - */ - bool PrunableSolver_(QueryTree *qnode, ReferenceTree *rnode, - Matrix ¤t_lanczos_vectors, - DRange &root_negative_dot_product_range, - DRange &root_positive_dot_product_range, - DRange &dsqd_range, - DRange &kernel_value_range, double &used_error); - - /** @brief The base-case exhaustive computation for dual-tree based - * computation of (B^T W(q) B) z(q). - * - * @param qnode The query node. - * @param rnode The reference node. - * @param current_lanczos_vectors Each column of this matrix is a current - * Lanczos vector for each query point. - */ - void DualtreeSolverBase_ - (QueryTree *qnode, ReferenceTree *rnode, const Matrix &qset, - const ArrayList &query_should_exit_the_loop, - const Matrix ¤t_lanczos_vectors, Matrix &lanczos_prod_l, - Matrix &lanczos_prod_e, Vector &lanczos_prod_used_error, - Vector &lanczos_prod_n_pruned, Matrix &neg_lanczos_prod_e, - Matrix &neg_lanczos_prod_u, Vector &neg_lanczos_prod_used_error, - Vector &neg_lanczos_prod_n_pruned); - - /** @brief The canonical case for dual-tree based computation of - * (B^T W(q) B) z(q) - * - * @param qnode The query node. - * @param rnode The reference node. - * @param current_lanczos_vectors Each column of this matrix is a current - * Lanczos vector for each query point. - */ - void DualtreeSolverCanonical_ - (QueryTree *qnode, ReferenceTree *rnode, const Matrix &qset, - const ArrayList &query_should_exit_the_loop, - const Matrix ¤t_lanczos_vectors, - Matrix &lanczos_prod_l, Matrix &lanczos_prod_e, - Vector &lanczos_prod_used_error, Vector &lanczos_prod_n_pruned, - Matrix &neg_lanczos_prod_e, Matrix &neg_lanczos_prod_u, - Vector &neg_lanczos_prod_used_error, Vector &neg_lanczos_prod_n_pruned); - - void SolveLeastSquaresByKrylov_(QueryTree *qroot, const Matrix &qset, - const Matrix &right_hand_sides, - Matrix &solution_vectors_e); - /** @brief Finalize the regression estimate for each query point by * taking the dot-product between [1; q^T] and the final * solution vector for (B^T W(q) B)^+ (B^T W(q) Y). @@ -471,13 +382,15 @@ class KrylovLpr { // point. This essentially becomes the right-hand side for each // query point. printf("Starting Phase 1...\n"); - ComputeRightHandSides_ - (qroot, qset, right_hand_sides_l, right_hand_sides_e, + ComputeWeightedVectorSum_ + (qroot, qset, rset_target_divided_by_norm_consts_, + right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); // TestRightHandSideComputation_(qset, right_hand_sides_e); printf("Phase 1 completed...\n"); + /* // The second phase solves the least squares problem: (B^T W(q) B) // z(q) = B^T W(q) Y for each query point q. printf("Starting Phase 2...\n"); @@ -504,6 +417,7 @@ class KrylovLpr { solution_vectors_e_single_alias); delete qroot_single; } + */ /* SolveLeastSquaresByKrylov_(qroot, qset, right_hand_sides_e, solution_vectors_e); @@ -568,7 +482,53 @@ class KrylovLpr { query_magnitude_weight_diagrams, (query_influence_values != NULL)); } - + + /** @brief Initialize the bandwidth by either fixed bandwidth + * parameter or a nearest neighbor based one (i.e. perform + * nearest neighbor and set the bandwidth equal to the k-th + * nearest neighbor distance). + */ + void InitializeBandwidths_() { + + kernels_.Init(rset_.n_cols()); + + if(fx_param_exists(NULL, "bandwidth")) { + printf("Using the fixed bandwidth method...\n"); + + double bandwidth = fx_param_double_req(NULL, "bandwidth"); + for(index_t i = 0; i < kernels_.size(); i++) { + kernels_[i].Init(bandwidth); + } + } + else { + printf("Using the nearest neighbor method...\n"); + AllkNN all_knn; + double knn_factor = fx_param_double(module_, "knn_factor", 0.2); + int knns = (int) (knn_factor * rset_.n_cols()); + all_knn.Init(rset_, 20, knns); + ArrayList resulting_neighbors; + ArrayList distances; + + all_knn.ComputeNeighbors(&resulting_neighbors, &distances); + + for(index_t i = 0; i < distances.size(); i += knns) { + kernels_[i / knns].Init(sqrt(distances[i + knns - 1])); + } + } + } + + void PrecomputeWeights_() { + rset_target_divided_by_norm_consts_.Init(rset_.n_cols()); + rset_inv_norm_consts_.Init(rset_.n_cols()); + + for(index_t i = 0; i < rset_.n_cols(); i++) { + rset_target_divided_by_norm_consts_[i] = + rset_targets_[i] / kernels_[i].CalcNormConstant(dimension_); + rset_inv_norm_consts_[i] = 1.0 / + kernels_[i].CalcNormConstant(dimension_); + } + } + public: ////////// Constructor/Destructor ////////// @@ -659,12 +619,9 @@ class KrylovLpr { rset_targets_.CopyValues(tmp_rset_targets); fx_timer_stop(NULL, "krylov_lpr_reference_tree_construct"); - // Initialize the kernel. - double bandwidth = fx_param_double_req(NULL, "bandwidth"); - kernels_.Init(rset_.n_cols()); - for(index_t i = 0; i < rset_.n_cols(); i++) { - kernels_[i].Init(bandwidth); - } + // Initialize the kernels. + InitializeBandwidths_(); + PrecomputeWeights_(); // Train the model using the reference set (i.e. compute // confidence interval and degrees of freedom.) @@ -672,7 +629,6 @@ class KrylovLpr { // initialize the reference side statistics. target_weighted_rset_.Init(row_length_, rset_.n_cols()); - InitializeReferenceStatistics_(rroot_); ComputeMain_(references, &rset_regression_estimates_, &rset_confidence_bands_, &rset_magnitude_weight_diagrams_, @@ -701,7 +657,6 @@ class KrylovLpr { }; #include "krylov_lpr_setup_impl.h" -#include "krylov_lpr_solver_impl.h" #include "krylov_lpr_test.h" #undef INSIDE_KRYLOV_LPR_H diff --git a/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h b/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h index d177145279..58c891fb4e 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h +++ b/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h @@ -9,7 +9,7 @@ template void KrylovLpr:: -InitializeQueryTreeRightHandSides_(QueryTree *qnode) { +InitializeQueryTree_(QueryTree *qnode) { // Set the bounds to default values. qnode->stat().Reset(); @@ -17,14 +17,15 @@ InitializeQueryTreeRightHandSides_(QueryTree *qnode) { // If the query node is not a leaf, then recurse. if(!qnode->is_leaf()) { - InitializeQueryTreeRightHandSides_(qnode->left()); - InitializeQueryTreeRightHandSides_(qnode->right()); + InitializeQueryTree_(qnode->left()); + InitializeQueryTree_(qnode->right()); } } template void KrylovLpr:: -InitializeReferenceStatistics_(ReferenceTree *rnode) { +InitializeReferenceStatistics_(ReferenceTree *rnode, int column_index, + const Vector &weights) { if(rnode->is_leaf()) { @@ -51,15 +52,11 @@ InitializeReferenceStatistics_(ReferenceTree *rnode) { // Compute the multiindex expansion of the given reference point. MultiIndexUtil::ComputePointMultivariatePolynomial (dimension_, lpr_order_, r_col.ptr(), reference_point_expansion.ptr()); - - // Accumulate each expansion onto its bounding box. - rnode->stat().reference_point_expansion_bound_ |= - reference_point_expansion; // Scale the expansion by the reference target. - la::ScaleOverwrite(row_length_, rset_targets_[r], - reference_point_expansion.ptr(), - r_target_weighted_by_coordinates); + la::ScaleOverwrite + (row_length_, weights[r] * reference_point_expansion[column_index], + reference_point_expansion.ptr(), r_target_weighted_by_coordinates); // Accumulate the far field coefficient for the target weighted // reference vector and the outerproduct. @@ -67,12 +64,6 @@ InitializeReferenceStatistics_(ReferenceTree *rnode) { rnode->stat().target_weighted_data_far_field_expansion_[j]. Add(r_target_weighted_by_coordinates[j], kernels_[r].bandwidth_sq(), r_col); - - for(index_t i = 0; i <= j; i++) { - rnode->stat().data_outer_products_far_field_expansion_[j][i]. - Add(reference_point_expansion[j] * reference_point_expansion[i], - kernels_[r].bandwidth_sq(), r_col); - } } // Tally up the weighted targets. @@ -97,8 +88,8 @@ InitializeReferenceStatistics_(ReferenceTree *rnode) { else { // Recursively call the function with left and right and merge. - InitializeReferenceStatistics_(rnode->left()); - InitializeReferenceStatistics_(rnode->right()); + InitializeReferenceStatistics_(rnode->left(), column_index, weights); + InitializeReferenceStatistics_(rnode->right(), column_index, weights); // Compute the sum of the sub sums. la::AddOverwrite((rnode->left()->stat()).sum_target_weighted_data_, @@ -117,26 +108,7 @@ InitializeReferenceStatistics_(ReferenceTree *rnode) { rnode->stat().target_weighted_data_far_field_expansion_[j]. Add(rnode->right()->stat(). target_weighted_data_far_field_expansion_[j]); - - for(index_t i = 0; i <= j; i++) { - - // First the far field moments of outer product using the bandwidth - rnode->stat().data_outer_products_far_field_expansion_[j][i]. - Add(rnode->left()->stat(). - data_outer_products_far_field_expansion_[j][i]); - rnode->stat().data_outer_products_far_field_expansion_[j][i]. - Add(rnode->right()->stat(). - data_outer_products_far_field_expansion_[j][i]); - } // end of iterating over each row. } // end of iterating over each column. - - // Combine the bounds of the reference point expansion owned by - // the two children. - rnode->stat().reference_point_expansion_bound_.Reset(); - rnode->stat().reference_point_expansion_bound_ |= - rnode->left()->stat().reference_point_expansion_bound_; - rnode->stat().reference_point_expansion_bound_ |= - rnode->right()->stat().reference_point_expansion_bound_; // Compute the min of the min bandwidths and the max of the max // bandwidths owned among the children. @@ -152,7 +124,7 @@ InitializeReferenceStatistics_(ReferenceTree *rnode) { } template -void KrylovLpr::DualtreeRightHandSidesBase_ +void KrylovLpr::DualtreeWeightedVectorSumBase_ (QueryTree *qnode, ReferenceTree *rnode, const Matrix &qset, Matrix &right_hand_sides_l, Matrix &right_hand_sides_e, Vector &right_hand_sides_used_error, Vector &right_hand_sides_n_pruned) { @@ -229,7 +201,7 @@ void KrylovLpr::DualtreeRightHandSidesBase_ } template -void KrylovLpr::DualtreeRightHandSidesCanonical_ +void KrylovLpr::DualtreeWeightedVectorSumCanonical_ (QueryTree *qnode, ReferenceTree *rnode, const Matrix &qset, Matrix &right_hand_sides_l, Matrix &right_hand_sides_e, Vector &right_hand_sides_used_error, Vector &right_hand_sides_n_pruned) { @@ -249,7 +221,7 @@ void KrylovLpr::DualtreeRightHandSidesCanonical_ (qnode, rnode, dsqd_range, kernel_value_range); // try finite difference pruning first - if(TPruneRule::PrunableKrylovRightHandSides + if(TPruneRule::PrunableWeightedVectorSum (internal_relative_error_, rnode->stat().sum_target_weighted_data_alloc_norm_, qnode, rnode, dsqd_range, kernel_value_range, delta_l, delta_e, @@ -290,7 +262,7 @@ void KrylovLpr::DualtreeRightHandSidesCanonical_ // for leaf pairs, go exhaustive if(rnode->is_leaf()) { - DualtreeRightHandSidesBase_ + DualtreeWeightedVectorSumBase_ (qnode, rnode, qset, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); return; @@ -301,10 +273,10 @@ void KrylovLpr::DualtreeRightHandSidesCanonical_ ReferenceTree *rnode_first = NULL, *rnode_second = NULL; LprUtil::BestReferenceNodePartners(qnode, rnode->left(), rnode->right(), &rnode_first, &rnode_second); - DualtreeRightHandSidesCanonical_ + DualtreeWeightedVectorSumCanonical_ (qnode, rnode_first, qset, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); - DualtreeRightHandSidesCanonical_ + DualtreeWeightedVectorSumCanonical_ (qnode, rnode_second, qset, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); return; @@ -345,10 +317,10 @@ void KrylovLpr::DualtreeRightHandSidesCanonical_ LprUtil::BestQueryNodePartners(rnode, qnode->left(), qnode->right(), &qnode_first, &qnode_second); - DualtreeRightHandSidesCanonical_ + DualtreeWeightedVectorSumCanonical_ (qnode_first, rnode, qset, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); - DualtreeRightHandSidesCanonical_ + DualtreeWeightedVectorSumCanonical_ (qnode_second, rnode, qset, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); } @@ -360,11 +332,11 @@ void KrylovLpr::DualtreeRightHandSidesCanonical_ LprUtil::BestReferenceNodePartners(qnode->left(), rnode->left(), rnode->right(), &rnode_first, &rnode_second); - DualtreeRightHandSidesCanonical_ + DualtreeWeightedVectorSumCanonical_ (qnode->left(), rnode_first, qset, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); - DualtreeRightHandSidesCanonical_ + DualtreeWeightedVectorSumCanonical_ (qnode->left(), rnode_second, qset, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); @@ -372,11 +344,11 @@ void KrylovLpr::DualtreeRightHandSidesCanonical_ LprUtil::BestReferenceNodePartners(qnode->right(), rnode->left(), rnode->right(), &rnode_first, &rnode_second); - DualtreeRightHandSidesCanonical_ + DualtreeWeightedVectorSumCanonical_ (qnode->right(), rnode_first, qset, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); - DualtreeRightHandSidesCanonical_ + DualtreeWeightedVectorSumCanonical_ (qnode->right(), rnode_second, qset, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); @@ -405,7 +377,7 @@ void KrylovLpr::DualtreeRightHandSidesCanonical_ } template -void KrylovLpr::FinalizeQueryTreeRightHandSides_ +void KrylovLpr::FinalizeQueryTree_ (QueryTree *qnode, const Matrix &qset, Matrix &right_hand_sides_l, Matrix &right_hand_sides_e, Vector &right_hand_sides_used_error, Vector &right_hand_sides_n_pruned) { @@ -472,11 +444,11 @@ void KrylovLpr::FinalizeQueryTreeRightHandSides_ (q_stat.postponed_moment_ll_vector_e_[i]); } - FinalizeQueryTreeRightHandSides_ - (qnode->left(), qset, right_hand_sides_l, right_hand_sides_e, - right_hand_sides_used_error, right_hand_sides_n_pruned); - FinalizeQueryTreeRightHandSides_ - (qnode->right(), qset, right_hand_sides_l, right_hand_sides_e, - right_hand_sides_used_error, right_hand_sides_n_pruned); + FinalizeQueryTree_(qnode->left(), qset, right_hand_sides_l, + right_hand_sides_e, right_hand_sides_used_error, + right_hand_sides_n_pruned); + FinalizeQueryTree_(qnode->right(), qset, right_hand_sides_l, + right_hand_sides_e, right_hand_sides_used_error, + right_hand_sides_n_pruned); } } diff --git a/fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h b/fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h deleted file mode 100644 index 9fe45b4547..0000000000 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h +++ /dev/null @@ -1,857 +0,0 @@ -// Make sure this file is included only in local_linear_krylov.h. This -// is not a public header file! -#ifndef INSIDE_KRYLOV_LPR_H -#error "This file is not a public header file!" -#endif - -template -void KrylovLpr::DualtreeSolverBase_ -(QueryTree *qnode, ReferenceTree *rnode, const Matrix &qset, - const ArrayList &query_should_exit_the_loop, - const Matrix ¤t_lanczos_vectors, - Matrix &lanczos_prod_l, Matrix &lanczos_prod_e, - Vector &lanczos_prod_used_error, Vector &lanczos_prod_n_pruned, - Matrix &neg_lanczos_prod_e, Matrix &neg_lanczos_prod_u, - Vector &neg_lanczos_prod_used_error, Vector &neg_lanczos_prod_n_pruned) { - - // Temporary variable for storing multivariate expansion of a - // reference point. - Vector reference_point_expansion; - reference_point_expansion.Init(row_length_); - - // Clear the summary statistics of the current query node so that we - // can refine it to better bounds. - qnode->stat().ll_vector_norm_l_ = DBL_MAX; - qnode->stat().ll_vector_used_error_ = 0; - qnode->stat().ll_vector_n_pruned_ = DBL_MAX; - qnode->stat().neg_ll_vector_norm_l_ = DBL_MAX; - qnode->stat().neg_ll_vector_used_error_ = 0; - qnode->stat().neg_ll_vector_n_pruned_ = DBL_MAX; - - // for each query point - for(index_t q = qnode->begin(); q < qnode->end(); q++) { - - // This is potentially inefficient and could be solved by - // rebuilding the query tree everytime when a query point exists - // the Lanczos outer loop. - if(query_should_exit_the_loop[q]) { - continue; - } - - // get query point. - const double *q_col = qset.GetColumnPtr(q); - - // Get the query point's associated current Lanczos vector. - const double *q_lanczos_vector = current_lanczos_vectors.GetColumnPtr(q); - - // get the column vectors accumulating the sums to update. - double *q_lanczos_prod_l = lanczos_prod_l.GetColumnPtr(q); - double *q_lanczos_prod_e = lanczos_prod_e.GetColumnPtr(q); - double *q_neg_lanczos_prod_e = neg_lanczos_prod_e.GetColumnPtr(q); - double *q_neg_lanczos_prod_u = neg_lanczos_prod_u.GetColumnPtr(q); - - // Incorporate the postponed information. - la::AddTo(row_length_, (qnode->stat().postponed_ll_vector_l_).ptr(), - q_lanczos_prod_l); - lanczos_prod_used_error[q] += - qnode->stat().postponed_ll_vector_used_error_; - lanczos_prod_n_pruned[q] += qnode->stat().postponed_ll_vector_n_pruned_; - - la::AddTo(row_length_, (qnode->stat().postponed_neg_ll_vector_u_).ptr(), - q_neg_lanczos_prod_u); - neg_lanczos_prod_used_error[q] += - qnode->stat().postponed_neg_ll_vector_used_error_; - neg_lanczos_prod_n_pruned[q] += - qnode->stat().postponed_neg_ll_vector_n_pruned_; - - // for each reference point - for(index_t r = rnode->begin(); r < rnode->end(); r++) { - - // get reference point. - const double *r_col = rset_.GetColumnPtr(r); - - // Compute the reference point expansion. - MultiIndexUtil::ComputePointMultivariatePolynomial - (dimension_, lpr_order_, r_col, reference_point_expansion.ptr()); - - // compute the pairwise squared distance and kernel value. - double dsqd = la::DistanceSqEuclidean(dimension_, q_col, r_col); - double kernel_value = kernels_[r].EvalUnnormOnSq(dsqd); - - // Take the dot product between the query point's Lanczos vector - // and the reference point expansion. - double dot_product = la::Dot(row_length_, q_lanczos_vector, - reference_point_expansion.ptr()); - double front_factor = dot_product * kernel_value; - - // For each vector component, update the lower/estimate/upper - // bound quantities. - if(front_factor > 0) { - la::AddExpert(row_length_, front_factor, - reference_point_expansion.ptr(), q_lanczos_prod_l); - la::AddExpert(row_length_, front_factor, - reference_point_expansion.ptr(), q_lanczos_prod_e); - } - else { - la::AddExpert(row_length_, front_factor, - reference_point_expansion.ptr(), q_neg_lanczos_prod_e); - la::AddExpert(row_length_, front_factor, - reference_point_expansion.ptr(), q_neg_lanczos_prod_u); - } - - } // end of iterating over each reference point. - - // Update the pruned quantities. - lanczos_prod_n_pruned[q] += - rnode->stat().sum_reference_point_expansion_norm_; - neg_lanczos_prod_n_pruned[q] += - rnode->stat().sum_reference_point_expansion_norm_; - - // Now, loop over each vector component for the current query and - // correct the upper bound by the assumption made in the - // initialization phase of the query tree. Refine min and max - // summary statistics. - qnode->stat().ll_vector_norm_l_ = - std::min(qnode->stat().ll_vector_norm_l_, - MatrixUtil::EntrywiseLpNorm(row_length_, q_lanczos_prod_l, 1)); - qnode->stat().ll_vector_used_error_ = - std::max(qnode->stat().ll_vector_used_error_, - lanczos_prod_used_error[q]); - qnode->stat().ll_vector_n_pruned_ = - std::min(qnode->stat().ll_vector_n_pruned_, lanczos_prod_n_pruned[q]); - - qnode->stat().neg_ll_vector_norm_l_ = - std::min(qnode->stat().neg_ll_vector_norm_l_, - MatrixUtil::EntrywiseLpNorm(row_length_, - q_neg_lanczos_prod_u, 1)); - qnode->stat().neg_ll_vector_used_error_ = - std::max(qnode->stat().neg_ll_vector_used_error_, - neg_lanczos_prod_used_error[q]); - qnode->stat().neg_ll_vector_n_pruned_ = - std::min(qnode->stat().neg_ll_vector_n_pruned_, - neg_lanczos_prod_n_pruned[q]); - - } // end of iterating over each query point. - - // Clear postponed information. - (qnode->stat().postponed_ll_vector_l_).SetZero(); - qnode->stat().postponed_ll_vector_used_error_ = 0; - qnode->stat().postponed_ll_vector_n_pruned_ = 0; - (qnode->stat().postponed_neg_ll_vector_u_).SetZero(); - qnode->stat().postponed_neg_ll_vector_used_error_ = 0; - qnode->stat().postponed_neg_ll_vector_n_pruned_ = 0; -} - -template -void KrylovLpr::DualtreeSolverCanonical_ -(QueryTree *qnode, ReferenceTree *rnode, const Matrix &qset, - const ArrayList &query_should_exit_the_loop, - const Matrix ¤t_lanczos_vectors, - Matrix &lanczos_prod_l, Matrix &lanczos_prod_e, - Vector &lanczos_prod_used_error, Vector &lanczos_prod_n_pruned, - Matrix &neg_lanczos_prod_e, Matrix &neg_lanczos_prod_u, - Vector &neg_lanczos_prod_used_error, Vector &neg_lanczos_prod_n_pruned) { - - // Variables for storing changes due to a prune. - double delta_used_error = 0, delta_n_pruned = 0, delta_neg_used_error = 0, - delta_neg_n_pruned = 0; - Vector delta_l, delta_e, delta_neg_u, delta_neg_e; - delta_l.Init(row_length_); - delta_e.Init(row_length_); - delta_neg_u.Init(row_length_); - delta_neg_e.Init(row_length_); - - // temporary variable for holding distance/kernel value bounds - DRange dsqd_range; - DRange kernel_value_range; - - // First compute distance/kernel value bounds and dot product bound - // ranges. - LprUtil::SqdistAndKernelRanges_(qnode, rnode, dsqd_range, - kernel_value_range); - - // Temporary variables hold the dot product ranges. - DRange negative_dot_product_range, positive_dot_product_range; - - // Compute the dot product range. - DotProductBetweenTwoBounds_(qnode, rnode, negative_dot_product_range, - positive_dot_product_range); - - // try finite difference pruning first - if(TPruneRule::PrunableKrylovSolver - (internal_relative_error_, - rnode->stat().sum_reference_point_expansion_norm_, - qnode, rnode, dsqd_range, kernel_value_range, - negative_dot_product_range, positive_dot_product_range, - delta_l, delta_e, delta_used_error, delta_n_pruned, - delta_neg_u, delta_neg_e, delta_neg_used_error, delta_neg_n_pruned)) { - - la::AddTo(delta_l, &(qnode->stat().postponed_ll_vector_l_)); - la::AddTo(delta_e, &(qnode->stat().postponed_ll_vector_e_)); - qnode->stat().postponed_ll_vector_used_error_ += delta_used_error; - qnode->stat().postponed_ll_vector_n_pruned_ += delta_n_pruned; - - la::AddTo(delta_neg_u, &(qnode->stat().postponed_neg_ll_vector_u_)); - la::AddTo(delta_neg_e, &(qnode->stat().postponed_neg_ll_vector_e_)); - qnode->stat().postponed_neg_ll_vector_used_error_ += delta_neg_used_error; - qnode->stat().postponed_neg_ll_vector_n_pruned_ += delta_neg_n_pruned; - - num_finite_difference_prunes_++; - - return; - } - - // For the Epanechnikov kernel, we can prune using the far field - // moments if the maximum distance between the two nodes is within - // the bandwidth! This if-statement does not apply to the Gaussian - // kernel, so I need to fix in the future! - if(rnode->stat().min_bandwidth_kernel.bandwidth_sq() >= dsqd_range.hi && - rnode->count() > dimension_ * dimension_) { - - la::AddTo(delta_l, &(qnode->stat().postponed_ll_vector_l_)); - qnode->stat().postponed_ll_vector_n_pruned_ += delta_n_pruned; - - la::AddTo(delta_neg_u, &(qnode->stat().postponed_neg_ll_vector_u_)); - qnode->stat().postponed_neg_ll_vector_n_pruned_ += delta_neg_n_pruned; - - // Add the Epanechnikov moments - for(index_t j = 0; j < row_length_; j++) { - for(index_t i = 0; i <= j; i++) { - qnode->stat().postponed_epanechnikov_moments_[j][i]. - Add(rnode->stat().data_outer_products_far_field_expansion_[j][i]); - } - } - - // Keep track of the far-field prunes. - num_epanechnikov_prunes_++; - return; - } - - // for leaf query node - if(qnode->is_leaf()) { - - // for leaf pairs, go exhaustive - if(rnode->is_leaf()) { - DualtreeSolverBase_ - (qnode, rnode, qset, query_should_exit_the_loop, - current_lanczos_vectors, lanczos_prod_l, lanczos_prod_e, - lanczos_prod_used_error, lanczos_prod_n_pruned, - neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - return; - } - - // for non-leaf reference, expand reference node - else { - ReferenceTree *rnode_first = NULL, *rnode_second = NULL; - LprUtil::BestReferenceNodePartners(qnode, rnode->left(), rnode->right(), - &rnode_first, &rnode_second); - DualtreeSolverCanonical_ - (qnode, rnode_first, qset, query_should_exit_the_loop, - current_lanczos_vectors, lanczos_prod_l, lanczos_prod_e, - lanczos_prod_used_error, lanczos_prod_n_pruned, - neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - DualtreeSolverCanonical_ - (qnode, rnode_second, qset, query_should_exit_the_loop, - current_lanczos_vectors, lanczos_prod_l, lanczos_prod_e, - lanczos_prod_used_error, lanczos_prod_n_pruned, - neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - - return; - } - } - - // for non-leaf query node - else { - - // Declare references to the query stats. - KrylovLprQStat &q_stat = qnode->stat(); - KrylovLprQStat &q_left_stat = qnode->left()->stat(); - KrylovLprQStat &q_right_stat = qnode->right()->stat(); - - // Push down postponed bound changes owned by the current query - // node to the children of the query node. - la::AddTo(q_stat.postponed_ll_vector_l_, - &(q_left_stat.postponed_ll_vector_l_)); - la::AddTo(q_stat.postponed_ll_vector_l_, - &(q_right_stat.postponed_ll_vector_l_)); - q_left_stat.postponed_ll_vector_used_error_ += - q_stat.postponed_ll_vector_used_error_; - q_right_stat.postponed_ll_vector_used_error_ += - q_stat.postponed_ll_vector_used_error_; - q_left_stat.postponed_ll_vector_n_pruned_ += - q_stat.postponed_ll_vector_n_pruned_; - q_right_stat.postponed_ll_vector_n_pruned_ += - q_stat.postponed_ll_vector_n_pruned_; - - la::AddTo(q_stat.postponed_neg_ll_vector_u_, - &(q_left_stat.postponed_neg_ll_vector_u_)); - la::AddTo(q_stat.postponed_neg_ll_vector_u_, - &(q_right_stat.postponed_neg_ll_vector_u_)); - q_left_stat.postponed_neg_ll_vector_used_error_ += - q_stat.postponed_neg_ll_vector_used_error_; - q_right_stat.postponed_neg_ll_vector_used_error_ += - q_stat.postponed_neg_ll_vector_used_error_; - q_left_stat.postponed_neg_ll_vector_n_pruned_ += - q_stat.postponed_neg_ll_vector_n_pruned_; - q_right_stat.postponed_neg_ll_vector_n_pruned_ += - q_stat.postponed_neg_ll_vector_n_pruned_; - - // Clear the statistics after pushing them downwards. - q_stat.postponed_ll_vector_l_.SetZero(); - q_stat.postponed_ll_vector_used_error_ = 0; - q_stat.postponed_ll_vector_n_pruned_ = 0; - q_stat.postponed_neg_ll_vector_u_.SetZero(); - q_stat.postponed_neg_ll_vector_used_error_ = 0; - q_stat.postponed_neg_ll_vector_n_pruned_ = 0; - - // For a leaf reference node, expand query node - if(rnode->is_leaf()) { - QueryTree *qnode_first = NULL, *qnode_second = NULL; - - LprUtil::BestQueryNodePartners(rnode, qnode->left(), qnode->right(), - &qnode_first, &qnode_second); - DualtreeSolverCanonical_ - (qnode_first, rnode, qset, query_should_exit_the_loop, - current_lanczos_vectors, lanczos_prod_l, lanczos_prod_e, - lanczos_prod_used_error, lanczos_prod_n_pruned, - neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - DualtreeSolverCanonical_ - (qnode_second, rnode, qset, query_should_exit_the_loop, - current_lanczos_vectors, lanczos_prod_l, lanczos_prod_e, - lanczos_prod_used_error, lanczos_prod_n_pruned, - neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - } - - // for non-leaf reference node, expand both query and reference nodes - else { - ReferenceTree *rnode_first = NULL, *rnode_second = NULL; - - LprUtil::BestReferenceNodePartners(qnode->left(), rnode->left(), - rnode->right(), &rnode_first, - &rnode_second); - DualtreeSolverCanonical_ - (qnode->left(), rnode_first, qset, query_should_exit_the_loop, - current_lanczos_vectors, lanczos_prod_l, lanczos_prod_e, - lanczos_prod_used_error, lanczos_prod_n_pruned, - neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - DualtreeSolverCanonical_ - (qnode->left(), rnode_second, qset, query_should_exit_the_loop, - current_lanczos_vectors, lanczos_prod_l, lanczos_prod_e, - lanczos_prod_used_error, lanczos_prod_n_pruned, - neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - - LprUtil::BestReferenceNodePartners(qnode->right(), rnode->left(), - rnode->right(), &rnode_first, - &rnode_second); - DualtreeSolverCanonical_ - (qnode->right(), rnode_first, qset, query_should_exit_the_loop, - current_lanczos_vectors, lanczos_prod_l, lanczos_prod_e, - lanczos_prod_used_error, lanczos_prod_n_pruned, - neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - - DualtreeSolverCanonical_ - (qnode->right(), rnode_second, qset, query_should_exit_the_loop, - current_lanczos_vectors, lanczos_prod_l, lanczos_prod_e, - lanczos_prod_used_error, lanczos_prod_n_pruned, - neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - } - - // reaccumulate the summary statistics. - q_stat.ll_vector_norm_l_ = - std::min - (q_left_stat.ll_vector_norm_l_ + - MatrixUtil::EntrywiseLpNorm(q_left_stat.postponed_ll_vector_l_, 1), - q_right_stat.ll_vector_norm_l_ + - MatrixUtil::EntrywiseLpNorm(q_right_stat.postponed_ll_vector_l_, 1)); - q_stat.ll_vector_used_error_ = - std::max(q_left_stat.ll_vector_used_error_ + - q_left_stat.postponed_ll_vector_used_error_, - q_right_stat.ll_vector_used_error_ + - q_right_stat.postponed_ll_vector_used_error_); - q_stat.ll_vector_n_pruned_ = - std::min(q_left_stat.ll_vector_n_pruned_ + - q_left_stat.postponed_ll_vector_n_pruned_, - q_right_stat.ll_vector_n_pruned_ + - q_right_stat.postponed_ll_vector_n_pruned_); - - q_stat.neg_ll_vector_norm_l_ = - std::min(q_left_stat.neg_ll_vector_norm_l_ + - MatrixUtil::EntrywiseLpNorm - (q_left_stat.postponed_neg_ll_vector_u_, 1), - q_right_stat.neg_ll_vector_norm_l_ + - MatrixUtil::EntrywiseLpNorm - (q_right_stat.postponed_neg_ll_vector_u_, 1)); - q_stat.neg_ll_vector_used_error_ = - std::max(q_left_stat.neg_ll_vector_used_error_ + - q_left_stat.postponed_neg_ll_vector_used_error_, - q_right_stat.neg_ll_vector_used_error_ + - q_right_stat.postponed_neg_ll_vector_used_error_); - q_stat.neg_ll_vector_n_pruned_ = - std::min(q_left_stat.neg_ll_vector_n_pruned_ + - q_left_stat.postponed_neg_ll_vector_n_pruned_, - q_right_stat.neg_ll_vector_n_pruned_ + - q_right_stat.postponed_neg_ll_vector_n_pruned_); - return; - } // end of the case: non-leaf query node. - -} - -template -void KrylovLpr::DotProductBetweenTwoBounds_ -(QueryTree *qnode, ReferenceTree *rnode, DRange &negative_dot_product_range, - DRange &positive_dot_product_range) { - - DHrectBound<2> lanczos_vectors_bound = qnode->stat().lanczos_vectors_bound_; - - // Initialize the dot-product ranges. - negative_dot_product_range.lo = negative_dot_product_range.hi = 0; - positive_dot_product_range.lo = positive_dot_product_range.hi = 0; - - for(index_t d = 0; d < row_length_; d++) { - - const DRange &lanczos_directional_bound = lanczos_vectors_bound.get(d); - const DRange &reference_node_directional_bound = - rnode->stat().reference_point_expansion_bound_.get(d); - - if(lanczos_directional_bound.lo > 0) { - positive_dot_product_range.lo += lanczos_directional_bound.lo * - reference_node_directional_bound.lo; - positive_dot_product_range.hi += lanczos_directional_bound.hi * - reference_node_directional_bound.hi; - } - else if(lanczos_directional_bound.Contains(0)) { - positive_dot_product_range.hi += lanczos_directional_bound.hi * - reference_node_directional_bound.hi; - negative_dot_product_range.lo += lanczos_directional_bound.lo * - reference_node_directional_bound.hi; - } - else { - negative_dot_product_range.lo += lanczos_directional_bound.lo * - reference_node_directional_bound.hi; - negative_dot_product_range.hi += lanczos_directional_bound.hi * - reference_node_directional_bound.lo; - } - } // End of looping over each component... -} - -template -void KrylovLpr::InitializeQueryTreeLanczosVectorBound_ -(QueryTree *qnode, const Matrix &qset, - const ArrayList &exclude_query_flag, - const Matrix ¤t_lanczos_vectors) { - - // Set the bound quantities to default. - qnode->stat().Reset(); - - // If the query node is a leaf, then exhaustively iterate over and - // form bounding boxes of the current solution. - if(qnode->is_leaf()) { - qnode->bound().Reset(); - - for(index_t q = qnode->begin(); q < qnode->end(); q++) { - - // If the current query point is not to be included in the - // bounding box, then skip it. - if(exclude_query_flag[q]) { - continue; - } - - Vector query_vector; - Vector lanczos_vector; - current_lanczos_vectors.MakeColumnVector(q, &lanczos_vector); - qset.MakeColumnVector(q, &query_vector); - qnode->stat().lanczos_vectors_bound_ |= lanczos_vector; - qnode->bound() |= query_vector; - } - } - - // Otherwise, traverse the left and the right and combine the - // bounding boxes of the solutions for the two children. - else { - InitializeQueryTreeLanczosVectorBound_(qnode->left(), qset, - exclude_query_flag, - current_lanczos_vectors); - InitializeQueryTreeLanczosVectorBound_(qnode->right(), qset, - exclude_query_flag, - current_lanczos_vectors); - - // Reset the bounding box for the Lanczos vectors and reform it - // using the bounding boxes owned by the children. - qnode->stat().lanczos_vectors_bound_ |= - (qnode->left()->stat()).lanczos_vectors_bound_; - qnode->stat().lanczos_vectors_bound_ |= - (qnode->right()->stat()).lanczos_vectors_bound_; - - // Ditto for the bounding box for the query points. - qnode->bound().Reset(); - qnode->bound() |= qnode->left()->bound(); - qnode->bound() |= qnode->right()->bound(); - } -} - -template -void KrylovLpr::FinalizeQueryTreeLanczosMultiplier_ -(QueryTree *qnode, const Matrix &qset, - const ArrayList &exclude_query_flag, - const Matrix ¤t_lanczos_vectors, - Matrix &lanczos_prod_l, Matrix &lanczos_prod_e, - Vector &lanczos_prod_used_error, Vector &lanczos_prod_n_pruned, - Matrix &neg_lanczos_prod_e, Matrix &neg_lanczos_prod_u, - Vector &neg_lanczos_prod_used_error, Vector &neg_lanczos_prod_n_pruned) { - - KrylovLprQStat &q_stat = qnode->stat(); - - if(qnode->is_leaf()) { - - // The matrix to store the evaluated moments at each query point. - Matrix evaluated_moments; - evaluated_moments.Init(row_length_, row_length_); - Vector evaluated_moments_times_lanczos_vector; - evaluated_moments_times_lanczos_vector.Init(row_length_); - - // Iterate over each query point. - for(index_t q = qnode->begin(); q < qnode->end(); q++) { - - if(exclude_query_flag[q]) { - continue; - } - - // Get the current query point. - Vector q_col; - qset.MakeColumnVector(q, &q_col); - - // Get the pointer to the current lanczos vector owned by the - // current query point. - Vector q_current_lanczos_vector; - current_lanczos_vectors.MakeColumnVector(q, &q_current_lanczos_vector); - - // Get the column vectors accumulating the sums to update. - double *q_lanczos_prod_l = lanczos_prod_l.GetColumnPtr(q); - double *q_lanczos_prod_e = lanczos_prod_e.GetColumnPtr(q); - double *q_neg_lanczos_prod_e = neg_lanczos_prod_e.GetColumnPtr(q); - double *q_neg_lanczos_prod_u = neg_lanczos_prod_u.GetColumnPtr(q); - - // Incorporate the postponed information. - la::AddTo(row_length_, (q_stat.postponed_ll_vector_l_).ptr(), - q_lanczos_prod_l); - la::AddTo(row_length_, (q_stat.postponed_ll_vector_e_).ptr(), - q_lanczos_prod_e); - la::AddTo(row_length_, (q_stat.postponed_neg_ll_vector_e_).ptr(), - q_neg_lanczos_prod_e); - la::AddTo(row_length_, (q_stat.postponed_neg_ll_vector_u_).ptr(), - q_neg_lanczos_prod_u); - - // Evaluate the Epanechnikov moments. - for(index_t i = 0; i < row_length_; i++) { - for(index_t j = 0; j <= i; j++) { - evaluated_moments.set - (j, i, qnode->stat().postponed_epanechnikov_moments_[i][j]. - ComputeKernelSum(q_col)); - } - } - for(index_t i = 0; i < row_length_; i++) { - for(index_t j = i + 1; j < row_length_; j++) { - evaluated_moments.set(j, i, evaluated_moments.get(i, j)); - } - } - - // Now compute the product between the evaluated moments and the - // Lanczos vector owned by this query point. - la::MulOverwrite(evaluated_moments, q_current_lanczos_vector, - &evaluated_moments_times_lanczos_vector); - - // Now accumulate the sum depending on the negativity or the - // positivity of each component. - for(index_t i = 0; i < row_length_; i++) { - if(evaluated_moments_times_lanczos_vector[i] > 0) { - q_lanczos_prod_e[i] += evaluated_moments_times_lanczos_vector[i]; - } - else { - q_neg_lanczos_prod_e[i] += evaluated_moments_times_lanczos_vector[i]; - } - } - - } // end of iterating over each query point. - } - else { - - KrylovLprQStat &q_left_stat = qnode->left()->stat(); - KrylovLprQStat &q_right_stat = qnode->right()->stat(); - - // Push down approximations - la::AddTo(q_stat.postponed_ll_vector_l_, - &(q_left_stat.postponed_ll_vector_l_)); - la::AddTo(q_stat.postponed_ll_vector_l_, - &(q_right_stat.postponed_ll_vector_l_)); - la::AddTo(q_stat.postponed_ll_vector_e_, - &(q_left_stat.postponed_ll_vector_e_)); - la::AddTo(q_stat.postponed_ll_vector_e_, - &(q_right_stat.postponed_ll_vector_e_)); - - la::AddTo(q_stat.postponed_neg_ll_vector_e_, - &(q_left_stat.postponed_neg_ll_vector_e_)); - la::AddTo(q_stat.postponed_neg_ll_vector_e_, - &(q_right_stat.postponed_neg_ll_vector_e_)); - la::AddTo(q_stat.postponed_neg_ll_vector_u_, - &(q_left_stat.postponed_neg_ll_vector_u_)); - la::AddTo(q_stat.postponed_neg_ll_vector_u_, - &(q_right_stat.postponed_neg_ll_vector_u_)); - - // Push down Epanechnikov pruned portions. - for(index_t i = 0; i < row_length_; i++) { - for(index_t j = 0; j <= i; j++) { - q_left_stat.postponed_epanechnikov_moments_[i][j].Add - (q_stat.postponed_epanechnikov_moments_[i][j]); - q_right_stat.postponed_epanechnikov_moments_[i][j].Add - (q_stat.postponed_epanechnikov_moments_[i][j]); - } - } - - // Recurse both branches of the query node. - FinalizeQueryTreeLanczosMultiplier_ - (qnode->left(), qset, exclude_query_flag, current_lanczos_vectors, - lanczos_prod_l, lanczos_prod_e, lanczos_prod_used_error, - lanczos_prod_n_pruned, neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - FinalizeQueryTreeLanczosMultiplier_ - (qnode->right(), qset, exclude_query_flag, current_lanczos_vectors, - lanczos_prod_l, lanczos_prod_e, lanczos_prod_used_error, - lanczos_prod_n_pruned, neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - } -} - -template -void KrylovLpr::SolveLeastSquaresByKrylov_ -(QueryTree *qroot, const Matrix &qset, const Matrix &right_hand_sides, - Matrix &solution_vectors_e) { - - // Initialize the initial solutions to be zero vectors. - solution_vectors_e.SetZero(); - - // Temporary variables needed for SYMMLQ iteration... - Matrix previous_lanczos_vectors; - Matrix current_lanczos_vectors; - Matrix v_tilde_mat; - previous_lanczos_vectors.Init(row_length_, qset.n_cols()); - v_tilde_mat.Init(row_length_, qset.n_cols()); - - Vector g_double_tilde_vec; - Vector g_vec; - g_double_tilde_vec.Init(qset.n_cols()); - g_vec.Init(qset.n_cols()); - - current_lanczos_vectors.Init(row_length_, qset.n_cols()); - - // More temporary variables for SYMMLQ routine... - Vector c_vec, beta_vec, beta_tilde_vec, s_vec; - Matrix w_mat; - c_vec.Init(qset.n_cols()); - beta_vec.Init(qset.n_cols()); - beta_tilde_vec.Init(qset.n_cols()); - s_vec.Init(qset.n_cols()); - w_mat.Init(row_length_, qset.n_cols()); - - // Initialize before entering the main iteration... This - // initialization implicitly assumes that initial guess to the - // linear system is the zero vector. - current_lanczos_vectors.CopyValues(right_hand_sides); - NormalizeMatrixColumnVectors_(current_lanczos_vectors, g_double_tilde_vec); - beta_vec.SetZero(); - beta_tilde_vec.SetZero(); - c_vec.SetAll(-1); - s_vec.SetZero(); - previous_lanczos_vectors.SetZero(); - w_mat.CopyValues(current_lanczos_vectors); - g_vec.SetZero(); - - // Flag to tell whether each query stays in the Krylov loop or not. - ArrayList query_should_exit_the_loop; - query_should_exit_the_loop.Init(qset.n_cols()); - - // Set the boolean flags to false - for(index_t q = 0; q < qset.n_cols(); q++) { - query_should_exit_the_loop[q] = false; - } - - // Initialize variables necessary for the dual-tree computation... - Matrix lanczos_prod_l, lanczos_prod_e, neg_lanczos_prod_e, - neg_lanczos_prod_u; - Vector lanczos_prod_used_error, lanczos_prod_n_pruned, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned; - lanczos_prod_l.Init(row_length_, qset.n_cols()); - lanczos_prod_e.Init(row_length_, qset.n_cols()); - neg_lanczos_prod_e.Init(row_length_, qset.n_cols()); - neg_lanczos_prod_u.Init(row_length_, qset.n_cols()); - lanczos_prod_used_error.Init(qset.n_cols()); - lanczos_prod_n_pruned.Init(qset.n_cols()); - neg_lanczos_prod_used_error.Init(qset.n_cols()); - neg_lanczos_prod_n_pruned.Init(qset.n_cols()); - - // Main iteration of the SYMMLQ algorithm - repeat until - // "convergence"... - for(index_t num_iter = 0; num_iter < sqrt(row_length_); num_iter++) { - - // Determine how many queries are in the Krylov loop. - int num_queries_in_krylov_loop = 0; - for(index_t q = 0; q < qset.n_cols(); q++) { - if(!query_should_exit_the_loop[q]) { - num_queries_in_krylov_loop++; - } - } - if(num_queries_in_krylov_loop == 0) { - break; - } - - // Initialize the query tree Lanzcos vector bounds. - InitializeQueryTreeLanczosVectorBound_ - (qroot, qset, query_should_exit_the_loop, current_lanczos_vectors); - - // Initialize the query tree bound statistics. - lanczos_prod_l.SetZero(); - lanczos_prod_e.SetZero(); - lanczos_prod_used_error.SetZero(); - lanczos_prod_n_pruned.SetZero(); - neg_lanczos_prod_e.SetZero(); - neg_lanczos_prod_u.SetZero(); - neg_lanczos_prod_used_error.SetZero(); - neg_lanczos_prod_n_pruned.SetZero(); - - // Multiply the current lanczos vector with the linear operator. - DualtreeSolverCanonical_ - (qroot, rroot_, qset, query_should_exit_the_loop, - current_lanczos_vectors, lanczos_prod_l, lanczos_prod_e, - lanczos_prod_used_error, lanczos_prod_n_pruned, neg_lanczos_prod_e, - neg_lanczos_prod_u, neg_lanczos_prod_used_error, - neg_lanczos_prod_n_pruned); - FinalizeQueryTreeLanczosMultiplier_ - (qroot, qset, query_should_exit_the_loop, current_lanczos_vectors, - lanczos_prod_l, lanczos_prod_e, lanczos_prod_used_error, - lanczos_prod_n_pruned, neg_lanczos_prod_e, neg_lanczos_prod_u, - neg_lanczos_prod_used_error, neg_lanczos_prod_n_pruned); - - // Compute v_tilde_mat (the residue after applying the linear - // operator the current Lanczos vector). - la::AddOverwrite(lanczos_prod_e, neg_lanczos_prod_e, &v_tilde_mat); - - /* - printf("Positive matrix: %g\n", - MatrixUtil::EntrywiseLpNorm(lanczos_prod_e, 1)); - printf("Negative matrix: %g\n", - MatrixUtil::EntrywiseLpNorm(neg_lanczos_prod_e, 1)); - TestKrylovComputation_(qset, v_tilde_mat, current_lanczos_vectors, - query_should_exit_the_loop); - */ - - for(index_t q = 0; q < qset.n_cols(); q++) { - - // If the current query is not in the Krylov loop, skip it. - if(query_should_exit_the_loop[q]) { - continue; - } - - double *v_tilde_mat_column = v_tilde_mat.GetColumnPtr(q); - double *previous_lanczos_vector = - previous_lanczos_vectors.GetColumnPtr(q); - double *current_lanczos_vector = - current_lanczos_vectors.GetColumnPtr(q); - - la::AddExpert(row_length_, -beta_vec[q], previous_lanczos_vector, - v_tilde_mat_column); - - // Compute alpha (a dot product b etween the current Lanczos - // vector and v_tilde vector). - double alpha = la::Dot(row_length_, current_lanczos_vector, - v_tilde_mat_column); - - // Subtract the component of the current Lanczos vector (a form - // of Gram-Schmidt orthogonalization.) - la::AddExpert(row_length_, -alpha, current_lanczos_vector, - v_tilde_mat_column); - - // Compute the length of v_tilde_mat_column and store into - // beta_vec. - beta_vec[q] = la::LengthEuclidean(row_length_, v_tilde_mat_column); - - // Make a backup copy of the current Lanczos vector. - for(index_t i = 0; i < row_length_; i++) { - previous_lanczos_vector[i] = current_lanczos_vector[i]; - } - - // Set a new current Lanczos vector based on v_tilde_mat_column. - // A potential place to watch out for division by zero!! - if(beta_vec[q] > 0) { - la::ScaleOverwrite(row_length_, 1.0 / beta_vec[q], v_tilde_mat_column, - current_lanczos_vector); - } - else { - query_should_exit_the_loop[q] = true; - la::ScaleOverwrite(row_length_, 1.0, v_tilde_mat_column, - current_lanczos_vector); - } - - // Compute l_1 - double l_1 = s_vec[q] * alpha - c_vec[q] * beta_tilde_vec[q]; - - // Compute l_2 - double l_2 = s_vec[q] * beta_vec[q]; - - // Compute alpha_tilde - double alpha_tilde = -s_vec[q] * beta_tilde_vec[q] - c_vec[q] * alpha; - - // Compute beta_tilde - beta_tilde_vec[q] = c_vec[q] * beta_vec[q]; - - double l_0 = sqrt(alpha_tilde * alpha_tilde + beta_vec[q] * beta_vec[q]); - - // Another potential place to watch for division by zero!! - if(l_0 != 0) { - c_vec[q] = alpha_tilde / l_0; - s_vec[q] = beta_vec[q] / l_0; - } - else { - query_should_exit_the_loop[q] = true; - } - - double g_tilde = g_double_tilde_vec[q] - l_1 * g_vec[q]; - g_double_tilde_vec[q] = -l_2 * g_vec[q]; - - // Another potential place to watch for division by zero!! - if(l_0 != 0) { - g_vec[q] = g_tilde / l_0; - } - else { - query_should_exit_the_loop[q] = true; - } - - // Update solution... - { - la::AddExpert(row_length_, g_vec[q] * c_vec[q], w_mat.GetColumnPtr(q), - solution_vectors_e.GetColumnPtr(q)); - la::AddExpert(row_length_, g_vec[q] * s_vec[q], current_lanczos_vector, - solution_vectors_e.GetColumnPtr(q)); - - la::Scale(row_length_, s_vec[q], w_mat.GetColumnPtr(q)); - la::AddExpert(row_length_, -c_vec[q], current_lanczos_vector, - w_mat.GetColumnPtr(q)); - } - - // Another criterion for quitting the Krylov loop... - if(sqrt(g_tilde * g_tilde + g_double_tilde_vec[q] * - g_double_tilde_vec[q]) < 0.001) { - query_should_exit_the_loop[q] = true; - } - - } // end of iterating over each query point. - - } // end of an iteration of SYMMLQ -} diff --git a/fastlib2/contrib/dongryel/regression/krylov_lpr_test.h b/fastlib2/contrib/dongryel/regression/krylov_lpr_test.h index b3914c90b2..e5d081b99f 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr_test.h +++ b/fastlib2/contrib/dongryel/regression/krylov_lpr_test.h @@ -33,7 +33,7 @@ void KrylovLpr::TestRightHandSideComputation_ // compute the pairwise squared distance and kernel value. double dsqd = la::DistanceSqEuclidean(dimension_, q_col, r_col); - double kernel_value = kernels_[0].EvalUnnormOnSq(dsqd); + double kernel_value = kernels_[r].EvalUnnormOnSq(dsqd); // Add up the contribution of the reference point. la::AddExpert(row_length_, kernel_value, r_weights, @@ -51,73 +51,3 @@ void KrylovLpr::TestRightHandSideComputation_ printf("Maximum relative error: %g\n", max_relative_error); } - -template -void KrylovLpr::TestKrylovComputation_ -(const Matrix &qset, const Matrix &approximated, - const Matrix ¤t_lanczos_vectors, - const ArrayList &query_should_exit_the_loop) { - - double max_relative_error = 0; - Matrix exact_vector_e; - exact_vector_e.Init(approximated.n_rows(), approximated.n_cols()); - exact_vector_e.SetZero(); - - Vector reference_point_expansion; - reference_point_expansion.Init(row_length_); - - for(index_t q = 0; q < qset.n_cols(); q++) { - - // If the current query should not be computed, then skip it. - if(query_should_exit_the_loop[q]) { - continue; - } - - // get the column vector corresponding to the current query point. - const double *q_col = qset.GetColumnPtr(q); - - // get the column vector corresponding to the Lanczos vector owned - // by the current query point. - const double *q_lanczos_vector = current_lanczos_vectors.GetColumnPtr(q); - - // get the column vector accumulating the sum. - Vector exact_vector_e_column, approx_column; - exact_vector_e.MakeColumnVector(q, &exact_vector_e_column); - approximated.MakeColumnVector(q, &approx_column); - - for(index_t r = 0; r < rset_.n_cols(); r++) { - - // get the column vector corresponding to the current reference point. - const double *r_col = rset_.GetColumnPtr(r); - - // compute the pairwise squared distance and kernel value. - double dsqd = la::DistanceSqEuclidean(dimension_, q_col, r_col); - double kernel_value = kernels_[0].EvalUnnormOnSq(dsqd); - - // Compute the reference point expansion. - MultiIndexUtil::ComputePointMultivariatePolynomial - (dimension_, lpr_order_, r_col, reference_point_expansion.ptr()); - - // Take the dot product between the query point's Lanczos vector - // and [1 r^T]^T. - double dot_product = - la::Dot(row_length_, q_lanczos_vector, - reference_point_expansion.ptr()); - double front_factor = dot_product * kernel_value; - - // Add the contribution of the current reference point. - la::AddExpert(row_length_, front_factor, reference_point_expansion.ptr(), - exact_vector_e_column.ptr()); - - } // end of iterating over each reference point. - - double relative_error = - MatrixUtil::EntrywiseNormDifferenceRelative - (exact_vector_e_column, approx_column, 1); - - max_relative_error = std::max(max_relative_error, relative_error); - - } // end of iterating over each query point. - - printf("Maximum relative error: %g\n", max_relative_error); -} diff --git a/fastlib2/contrib/dongryel/regression/krylov_stat.h b/fastlib2/contrib/dongryel/regression/krylov_stat.h index f2e887bcdb..da2c408c1d 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_stat.h +++ b/fastlib2/contrib/dongryel/regression/krylov_stat.h @@ -35,22 +35,6 @@ class KrylovLprRStat { */ ArrayList< EpanKernelMomentInfo > target_weighted_data_far_field_expansion_; - - /** @brief The far field expansion created by the outer - * products. The (i, j)-th element denotes the - * far-field expansion of the (i, j)-th component of - * the sum_data_outer_products_ matrix. - */ - ArrayList< ArrayList< EpanKernelMomentInfo > > - data_outer_products_far_field_expansion_; - - /** @brief The vector summing up the reference point expansion. - */ - Vector sum_reference_point_expansion_; - - /** @brief The norm of the sum_reference_point_expansion_ - */ - double sum_reference_point_expansion_norm_; /** @brief The minimum bandwidth among the reference point. */ @@ -60,9 +44,6 @@ class KrylovLprRStat { */ TKernel max_bandwidth_kernel; - /** @brief The bounding box for the reference point expansion */ - DHrectBound<2> reference_point_expansion_bound_; - ////////// Constructor/Destructor ////////// /** @brief The constructor which does not do anything. */ @@ -81,15 +62,10 @@ class KrylovLprRStat { sum_target_weighted_data_error_norm_ = 0; sum_target_weighted_data_alloc_norm_ = 0; - sum_reference_point_expansion_.SetZero(); - sum_reference_point_expansion_norm_ = 0; - // Initialize the bandwidth information to defaults. min_bandwidth_kernel.Init(DBL_MAX); max_bandwidth_kernel.Init(0); - reference_point_expansion_bound_.Reset(); - for(index_t j = 0; j < target_weighted_data_far_field_expansion_.size(); j++) { target_weighted_data_far_field_expansion_[j].Reset(); @@ -110,21 +86,10 @@ class KrylovLprRStat { sum_target_weighted_data_.Init(matrix_dimension); target_weighted_data_far_field_expansion_.Init(matrix_dimension); - sum_reference_point_expansion_.Init(matrix_dimension); for(index_t j = 0; j < matrix_dimension; j++) { target_weighted_data_far_field_expansion_[j].Init(dimension); } - - // Initialize memory for bound on outer product expansions. - reference_point_expansion_bound_.Init(matrix_dimension); - data_outer_products_far_field_expansion_.Init(matrix_dimension); - for(index_t j = 0; j < matrix_dimension; j++) { - data_outer_products_far_field_expansion_[j].Init(j + 1); - for(index_t i = 0; i <= j; i++) { - data_outer_products_far_field_expansion_[j][i].Init(dimension); - } - } } /** @brief Computing the statistics for a leaf node involves @@ -169,22 +134,6 @@ public: */ double ll_vector_n_pruned_; - /** @brief The lower bound on the norm of the negative components - * of the vector computation. - */ - double neg_ll_vector_norm_l_; - - /** @brief The upper bound on the used error for approximating the - * negative components of the vector computation. - */ - double neg_ll_vector_used_error_; - - /** @brief The lower bound on the portion of the reference set - * pruned for the query points owned by this node for the - * negative components. - */ - double neg_ll_vector_n_pruned_; - /** @brief The lower bound vector offset passed from the above on * each sum component of the vector owned by this node. */ @@ -208,35 +157,6 @@ public: */ double postponed_ll_vector_n_pruned_; - /** @brief This stores the portion pruned by finite difference for - * each negative sum component of the vector owned by this - * node. - */ - Vector postponed_neg_ll_vector_e_; - - /** @brief The upper bound vector offset passed from above on each - * negative sum component of the right hand sides owned by - * this node. - */ - Vector postponed_neg_ll_vector_u_; - - /** @brief The amount of used error passed down from above for - * approximating the negative components of the vector sum. - */ - double postponed_neg_ll_vector_used_error_; - - /** @brief The portion of the reference set pruned for approximating - * the negative components of the vector sum passed down - * from above. - */ - double postponed_neg_ll_vector_n_pruned_; - - /** @brief The bounding box for the Lanczos vectors. */ - DHrectBound<2> lanczos_vectors_bound_; - - ArrayList< ArrayList < EpanKernelMomentInfo > > - postponed_epanechnikov_moments_; - ////////// Constructor/Destructor ////////// /** @brief The constructor which does not do anything. */ @@ -253,24 +173,13 @@ public: ll_vector_norm_l_ = 0; ll_vector_used_error_ = 0; ll_vector_n_pruned_ = 0; - neg_ll_vector_norm_l_ = 0; - neg_ll_vector_used_error_ = 0; - neg_ll_vector_n_pruned_ = 0; postponed_ll_vector_l_.SetZero(); postponed_ll_vector_e_.SetZero(); postponed_ll_vector_used_error_ = 0; postponed_ll_vector_n_pruned_ = 0; - postponed_neg_ll_vector_e_.SetZero(); - postponed_neg_ll_vector_u_.SetZero(); - postponed_neg_ll_vector_used_error_ = 0; - postponed_neg_ll_vector_n_pruned_ = 0; - lanczos_vectors_bound_.Reset(); for(index_t i = 0; i < postponed_moment_ll_vector_e_.size(); i++) { postponed_moment_ll_vector_e_[i].Reset(); - for(index_t j = 0; j <= i; j++) { - postponed_epanechnikov_moments_[i][j].Reset(); - } } } @@ -289,18 +198,9 @@ public: postponed_ll_vector_l_.Init(matrix_dimension); postponed_ll_vector_e_.Init(matrix_dimension); postponed_moment_ll_vector_e_.Init(matrix_dimension); - postponed_epanechnikov_moments_.Init(matrix_dimension); for(index_t i = 0; i < postponed_moment_ll_vector_e_.size(); i++) { postponed_moment_ll_vector_e_[i].Init(dimension); - postponed_epanechnikov_moments_[i].Init(i + 1); - for(index_t j = 0; j <= i; j++) { - postponed_epanechnikov_moments_[i][j].Init(dimension); - } } - postponed_neg_ll_vector_e_.Init(matrix_dimension); - postponed_neg_ll_vector_u_.Init(matrix_dimension); - - lanczos_vectors_bound_.Init(matrix_dimension); } /** @brief Computing the statistics for a leaf node involves diff --git a/fastlib2/contrib/dongryel/regression/relative_prune_lpr.h b/fastlib2/contrib/dongryel/regression/relative_prune_lpr.h index 04849245ac..4220573eb5 100644 --- a/fastlib2/contrib/dongryel/regression/relative_prune_lpr.h +++ b/fastlib2/contrib/dongryel/regression/relative_prune_lpr.h @@ -153,7 +153,7 @@ class RelativePruneLpr { } template - static bool PrunableKrylovRightHandSides + static bool PrunableWeightedVectorSum (double relative_error, double total_alloc_error, QueryTree *qnode, ReferenceTree *rnode, const DRange &dsqd_range, const DRange &kernel_value_range, Vector &delta_l, Vector &delta_e, @@ -192,92 +192,6 @@ class RelativePruneLpr { // check pruning condition return (delta_used_error <= allowed_err); } - - template - static bool PrunableKrylovSolver - (double relative_error, double total_alloc_error, QueryTree *qnode, - ReferenceTree *rnode, const DRange &dsqd_range, - const DRange &kernel_value_range, - const DRange &negative_dot_product_range, - const DRange &positive_dot_product_range, - Vector &delta_l, Vector &delta_e, double &delta_used_error, - double &delta_n_pruned, Vector &delta_neg_u, Vector &delta_neg_e, - double &delta_neg_used_error, double &delta_neg_n_pruned) { - - // Compute the vector component lower and upper bound changes. This - // assumes that the maximum kernel value is 1. - la::ScaleOverwrite(positive_dot_product_range.lo * kernel_value_range.lo, - rnode->stat().sum_reference_point_expansion_, - &delta_l); - la::ScaleOverwrite(0.5 * (positive_dot_product_range.lo * - kernel_value_range.lo + - positive_dot_product_range.hi * - kernel_value_range.hi), - rnode->stat().sum_reference_point_expansion_, - &delta_e); - - la::ScaleOverwrite(0.5 * (negative_dot_product_range.lo * - kernel_value_range.hi + - negative_dot_product_range.hi * - kernel_value_range.lo), - rnode->stat().sum_reference_point_expansion_, - &delta_neg_e); - la::ScaleOverwrite(negative_dot_product_range.hi * kernel_value_range.lo, - rnode->stat().sum_reference_point_expansion_, - &delta_neg_u); - - // Compute the L1 norm of the most refined lower bound. - double new_ll_vector_norm_l = - qnode->stat().ll_vector_norm_l_ + - MatrixUtil::EntrywiseLpNorm(qnode->stat().postponed_ll_vector_l_, 1) + - MatrixUtil::EntrywiseLpNorm(delta_l, 1); - double new_ll_vector_used_error = qnode->stat().ll_vector_used_error_ + - qnode->stat().postponed_ll_vector_used_error_; - double new_ll_vector_n_pruned = qnode->stat().ll_vector_n_pruned_ + - qnode->stat().postponed_ll_vector_n_pruned_; - - double new_neg_ll_vector_norm_l = - qnode->stat().neg_ll_vector_norm_l_ + - MatrixUtil::EntrywiseLpNorm - (qnode->stat().postponed_neg_ll_vector_u_, 1) + - MatrixUtil::EntrywiseLpNorm(delta_neg_u, 1); - double new_neg_ll_vector_used_error = - qnode->stat().neg_ll_vector_used_error_ + - qnode->stat().postponed_neg_ll_vector_used_error_; - double new_neg_ll_vector_n_pruned = - qnode->stat().neg_ll_vector_n_pruned_ + - qnode->stat().postponed_neg_ll_vector_n_pruned_; - - // Compute the allowed amount of error for pruning the given query - // and reference pair. - double allowed_err = - (relative_error * new_ll_vector_norm_l - new_ll_vector_used_error) / - (total_alloc_error - new_ll_vector_n_pruned); - double neg_allowed_err = - (relative_error * new_neg_ll_vector_norm_l - - new_neg_ll_vector_used_error) / - (total_alloc_error - new_neg_ll_vector_n_pruned); - - // Record how much error and pruned portion will be if pruning - // were to succeed. - delta_used_error = - 0.5 * (positive_dot_product_range.hi * kernel_value_range.hi - - positive_dot_product_range.lo * kernel_value_range.lo) * - (rnode->stat().sum_reference_point_expansion_norm_); - delta_n_pruned = - rnode->stat().sum_reference_point_expansion_norm_; - - delta_neg_used_error = - 0.5 * (negative_dot_product_range.hi * kernel_value_range.lo - - negative_dot_product_range.lo * kernel_value_range.hi) * - (rnode->stat().sum_reference_point_expansion_norm_); - delta_neg_n_pruned = - rnode->stat().sum_reference_point_expansion_norm_; - - // check pruning condition - return (delta_used_error <= allowed_err && - delta_neg_used_error <= neg_allowed_err); - } }; #endif