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

Removed the custom-built Krylov solver, now going to replace it with Epetra solver

Browse Source
This commit is contained in:
Dongryeol Lee
2008-02-24 21:52:16 +00:00
parent 2b5b4a2e0e
commit fb5adf58b7
7 changed files with 121 additions and 1308 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib2/contrib/dongryel/regression/build.py
-1
View File
@@ -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"],
fastlib2/contrib/dongryel/regression/krylov_lpr.h
+90 -135
View File
@@ -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 &current_lanczos_vectors,
const ArrayList<bool> &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<bool> &exclude_query_flag,
const Matrix &current_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<bool> &exclude_query_flag,
const Matrix &current_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 &current_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<bool> &query_should_exit_the_loop,
const Matrix &current_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<bool> &query_should_exit_the_loop,
const Matrix &current_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<index_t> resulting_neighbors;
ArrayList<double> 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
fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h
+29 -57
View File
@@ -9,7 +9,7 @@
template<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::
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<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::
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<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::DualtreeRightHandSidesBase_
void KrylovLpr<TKernel, TPruneRule>::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<TKernel, TPruneRule>::DualtreeRightHandSidesBase_
}
template<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::DualtreeRightHandSidesCanonical_
void KrylovLpr<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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<TKernel, TPruneRule>::DualtreeRightHandSidesCanonical_
}
template<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::FinalizeQueryTreeRightHandSides_
void KrylovLpr<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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);
}
}
fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h
-857
View File
@@ -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<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::DualtreeSolverBase_
(QueryTree *qnode, ReferenceTree *rnode, const Matrix &qset,
const ArrayList<bool> &query_should_exit_the_loop,
const Matrix &current_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<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::DualtreeSolverCanonical_
(QueryTree *qnode, ReferenceTree *rnode, const Matrix &qset,
const ArrayList<bool> &query_should_exit_the_loop,
const Matrix &current_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<TKernel> &q_stat = qnode->stat();
KrylovLprQStat<TKernel> &q_left_stat = qnode->left()->stat();
KrylovLprQStat<TKernel> &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<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::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<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::InitializeQueryTreeLanczosVectorBound_
(QueryTree *qnode, const Matrix &qset,
const ArrayList<bool> &exclude_query_flag,
const Matrix &current_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<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::FinalizeQueryTreeLanczosMultiplier_
(QueryTree *qnode, const Matrix &qset,
const ArrayList<bool> &exclude_query_flag,
const Matrix &current_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<TKernel> &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<TKernel> &q_left_stat = qnode->left()->stat();
KrylovLprQStat<TKernel> &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<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::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<bool> 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
}
fastlib2/contrib/dongryel/regression/krylov_lpr_test.h
+1 -71
View File
@@ -33,7 +33,7 @@ void KrylovLpr<TKernel, TPruneRule>::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<TKernel, TPruneRule>::TestRightHandSideComputation_
printf("Maximum relative error: %g\n", max_relative_error);
}
template<typename TKernel, typename TPruneRule>
void KrylovLpr<TKernel, TPruneRule>::TestKrylovComputation_
(const Matrix &qset, const Matrix &approximated,
const Matrix &current_lanczos_vectors,
const ArrayList<bool> &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);
}
fastlib2/contrib/dongryel/regression/krylov_stat.h
-100
View File
@@ -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
fastlib2/contrib/dongryel/regression/relative_prune_lpr.h
+1 -87
View File
@@ -153,7 +153,7 @@ class RelativePruneLpr {
}
template<typename QueryTree, typename ReferenceTree>
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<typename QueryTree, typename ReferenceTree>
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