diff --git a/fastlib2/contrib/dongryel/regression/krylov_lpr.h b/fastlib2/contrib/dongryel/regression/krylov_lpr.h index b120a2ccce..bf31b3358f 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr.h +++ b/fastlib2/contrib/dongryel/regression/krylov_lpr.h @@ -165,8 +165,11 @@ class KrylovLpr { * ComputeRightHandSides_ function call to test the * correctness. */ - void TestRightHandSideComputation_(const Matrix &qset, - const Matrix &approximated); + void TestDualtreeComputation_ + (const Matrix &qset, const ArrayList *query_in_cg_loop, + const bool confidence_band_computation_phase, + const Vector &reference_weights, index_t column_index, + const Matrix &approximated); /** @brief Initialize the bound statistics relevant to the right * hand side computation. @@ -319,16 +322,16 @@ class KrylovLpr { rset_regression_estimates_[i]; rset_variance_ += prediction_error * prediction_error; } - - // This could happen if enough matrices are singular... + + rset_variance_ *= 1.0 / + (rset_.n_cols() - 2.0 * rset_first_degree_of_freedom_ + + rset_second_degree_of_freedom_); + + // This MIGHT happen if we have too few data points... if(rset_.n_cols() - 2.0 * rset_first_degree_of_freedom_ + rset_second_degree_of_freedom_ <= 0) { rset_variance_ = DBL_MAX; } - - rset_variance_ *= 1.0 / - (rset_.n_cols() - 2.0 * rset_first_degree_of_freedom_ + - rset_second_degree_of_freedom_); } void ComputeConfidenceBands_(const Matrix &queries, @@ -402,7 +405,7 @@ class KrylovLpr { InitializeQueryTree_(qroot, qset, query_in_cg_loop); // Call dualtree function. - if(query_in_cg_loop == NULL) { + if(query_in_cg_loop == NULL && !confidence_band_computation_phase) { DualtreeWeightedVectorSumCanonical_ (qroot, rroot_, qset, query_in_cg_loop, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, @@ -481,7 +484,15 @@ class KrylovLpr { (qroot, qset, rset_target_divided_by_norm_consts_, NULL, false, 0, right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned, leave_one_out_right_hand_sides_e); - + + // Uncomment the following three lines to test the correctness of + // the weighted vector sum. + /* + TestDualtreeComputation_(qset, NULL, false, + rset_target_divided_by_norm_consts_, 0, + right_hand_sides_e); + */ + printf("Phase 1 completed...\n"); // The second phase solves the least squares problem: (B^T W(q) B) @@ -646,25 +657,14 @@ class KrylovLpr { rset_inv_norm_consts_.Init(rset_.n_cols()); rset_inv_squared_norm_consts_.Init(rset_.n_cols()); - // Find out the minimum normalization constant. This assumes that - // the minimum normalization constant is greater than zero... - double min_norm_const = DBL_MAX; - for(index_t i = 0; i < rset_.n_cols(); i++) { - min_norm_const = std::min(min_norm_const, - kernels_[i].CalcNormConstant(dimension_)); - } - min_norm_const = 1; - for(index_t i = 0; i < rset_.n_cols(); i++) { rset_target_divided_by_norm_consts_[i] = - rset_targets_[i] / - (kernels_[i].CalcNormConstant(dimension_) / min_norm_const); - rset_inv_norm_consts_[i] = 1.0 / - (kernels_[i].CalcNormConstant(dimension_) / min_norm_const); + rset_targets_[i] / kernels_[i].CalcNormConstant(dimension_); + rset_inv_norm_consts_[i] = + 1.0 / kernels_[i].CalcNormConstant(dimension_); rset_inv_squared_norm_consts_[i] = 1.0 / (kernels_[i].CalcNormConstant(dimension_) * - kernels_[i].CalcNormConstant(dimension_) / - (min_norm_const * min_norm_const)); + kernels_[i].CalcNormConstant(dimension_)); } } diff --git a/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h b/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h index c1ae2cb548..ae0aa8358d 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h +++ b/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h @@ -78,7 +78,7 @@ void KrylovLpr::InitializeReferenceStatistics_ 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. for(index_t j = 0; j < row_length_; j++) { @@ -205,7 +205,7 @@ void KrylovLpr::DualtreeWeightedVectorSumBase_ // quantities. la::AddExpert(row_length_, kernel_value, r_weights, q_right_hand_side_l); la::AddExpert(row_length_, kernel_value, r_weights, q_right_hand_side_e); - + } // end of iterating over each reference point. // The current query point now has taken care of all reference @@ -621,8 +621,8 @@ void KrylovLpr::DecideComputationMethod_ // Otherwise, if we cannot prune, then exhaustively compute. DualtreeWeightedVectorSumBase_ - (qnode, rnode, qset, query_in_cg_loop, true, right_hand_sides_l, - right_hand_sides_e, right_hand_sides_used_error, + (qnode, rnode, qset, query_in_cg_loop, confidence_band_computation_phase, + right_hand_sides_l, right_hand_sides_e, right_hand_sides_used_error, right_hand_sides_n_pruned); } @@ -669,7 +669,7 @@ void KrylovLpr::StratifiedComputation_ // If the current query point is not in the CG loop, then skip // it. - if(!((*query_in_cg_loop)[q])) { + if(query_in_cg_loop != NULL && !((*query_in_cg_loop)[q])) { continue; } diff --git a/fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h b/fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h index 620cff7318..31b4753089 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h +++ b/fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h @@ -30,6 +30,11 @@ void KrylovLpr::LinearOperatorConfidenceBand (qroot, qset, rset_inv_squared_norm_consts_, NULL, true, d, vector_l, vector_e, vector_used_error, vector_n_pruned, NULL); + /* + TestDualtreeComputation_(qset, NULL, true, + rset_inv_squared_norm_consts_, d, vector_e); + */ + // Accumulate the product between the computed vector and each // scalar component of the X. for(index_t q = 0; q < qset.n_cols(); q++) { @@ -78,6 +83,11 @@ void KrylovLpr::LinearOperator vector_l, vector_e, vector_used_error, vector_n_pruned, leave_one_out_vector_e); + /* + TestDualtreeComputation_(qset, &query_in_cg_loop, false, + rset_inv_norm_consts_, d, vector_e); + */ + // Accumulate the product between the computed vector and each // scalar component of the X. for(index_t q = 0; q < qset.n_cols(); q++) { diff --git a/fastlib2/contrib/dongryel/regression/krylov_lpr_test.h b/fastlib2/contrib/dongryel/regression/krylov_lpr_test.h index e5d081b99f..45685c5ba4 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr_test.h +++ b/fastlib2/contrib/dongryel/regression/krylov_lpr_test.h @@ -5,9 +5,16 @@ #endif template -void KrylovLpr::TestRightHandSideComputation_ -(const Matrix &qset, const Matrix &approximated) { - +void KrylovLpr::TestDualtreeComputation_ +(const Matrix &qset, const ArrayList *query_in_cg_loop, + const bool confidence_band_computation_phase, + const Vector &reference_weights, index_t column_index, + const Matrix &approximated) { + + // temporary space for storing reference point expansion + Vector r_col_expansion; + r_col_expansion.Init(row_length_); + Matrix exact_vector_e; exact_vector_e.Init(approximated.n_rows(), approximated.n_cols()); exact_vector_e.SetZero(); @@ -15,6 +22,10 @@ void KrylovLpr::TestRightHandSideComputation_ for(index_t q = 0; q < qset.n_cols(); q++) { + if(query_in_cg_loop != NULL && !((*query_in_cg_loop)[q])) { + continue; + } + // get the column vector corresponding to the current query point. const double *q_col = qset.GetColumnPtr(q); @@ -28,19 +39,25 @@ void KrylovLpr::TestRightHandSideComputation_ // get the column vector corresponding to the current reference point. const double *r_col = rset_.GetColumnPtr(r); - // get the column vector containing the appropriate weights. - const double *r_weights = target_weighted_rset_.GetColumnPtr(r); + // compute the reference point expansion + MultiIndexUtil::ComputePointMultivariatePolynomial + (dimension_, lpr_order_, r_col, r_col_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); + + if(confidence_band_computation_phase) { + kernel_value *= kernel_value; + } // Add up the contribution of the reference point. - la::AddExpert(row_length_, kernel_value, r_weights, + la::AddExpert(row_length_, kernel_value * reference_weights[r] * + r_col_expansion[column_index], r_col_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); diff --git a/fastlib2/contrib/dongryel/regression/naive_lpr.h b/fastlib2/contrib/dongryel/regression/naive_lpr.h index 6f4df72afa..3ee6b872e3 100644 --- a/fastlib2/contrib/dongryel/regression/naive_lpr.h +++ b/fastlib2/contrib/dongryel/regression/naive_lpr.h @@ -328,16 +328,16 @@ class NaiveLpr { rset_variance_ += prediction_error * prediction_error; } + rset_variance_ *= 1.0 / + (rset_.n_cols() - 2.0 * rset_first_degree_of_freedom_ + + rset_second_degree_of_freedom_); + // This MIGHT happen if we have too few data points... if(rset_.n_cols() - 2.0 * rset_first_degree_of_freedom_ + rset_second_degree_of_freedom_ <= 0) { rset_variance_ = DBL_MAX; } - rset_variance_ *= 1.0 / - (rset_.n_cols() - 2.0 * rset_first_degree_of_freedom_ + - rset_second_degree_of_freedom_); - fx_format_result(module_, "reference_set_first_degree_of_freedom", "%g", rset_first_degree_of_freedom_); fx_format_result(module_, "reference_set_second_degree_of_freedom", @@ -474,7 +474,8 @@ class NaiveLpr { fx_timer_start(module_, "naive_lpr_querying_time"); ComputeMain_(queries, query_regression_estimates, NULL, - query_confidence_bands, query_magnitude_weight_diagrams); + query_confidence_bands, query_magnitude_weight_diagrams, + NULL); fx_timer_stop(module_, "naive_lpr_querying_time"); }