From f6dba35e8b73a1e81ca58b2a8884fceaa183e46d Mon Sep 17 00:00:00 2001 From: Dongryeol Lee Date: Sun, 2 Mar 2008 04:39:00 +0000 Subject: [PATCH] Fixed a bug in the dense version of the algorithm which did not shuffle the computed magnitude of weight diagrams and influence values accordingly due to the query tree construction; also there was a logic mistake in the design of the code which assumed that under variable-bandwidth case, the epanechnikov moment computed for the denominator matrix could be used for the weight diagram matrix, but this is not true since the effective normalization constant for weight diagram matrix requires squared normalization constants --- .../contrib/dongryel/regression/dense_lpr.h | 23 ++++- .../dongryel/regression/dense_lpr_impl.h | 85 +++++++++++++------ .../contrib/dongryel/regression/krylov_lpr.h | 10 ++- .../regression/krylov_lpr_setup_impl.h | 2 +- .../regression/krylov_lpr_solver_impl.h | 2 +- .../contrib/dongryel/regression/naive_lpr.h | 8 +- 6 files changed, 97 insertions(+), 33 deletions(-) diff --git a/fastlib2/contrib/dongryel/regression/dense_lpr.h b/fastlib2/contrib/dongryel/regression/dense_lpr.h index 0d57493a23..bdce544f42 100644 --- a/fastlib2/contrib/dongryel/regression/dense_lpr.h +++ b/fastlib2/contrib/dongryel/regression/dense_lpr.h @@ -75,6 +75,14 @@ class DenseLpr { ArrayList< ArrayList< EpanKernelMomentInfo > > data_outer_products_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 > > + weight_diagram_far_field_expansion_; + /** @brief The vector summing up the reference polynomial term * weighted by its target training value (i.e. B^T Y). */ @@ -118,6 +126,7 @@ class DenseLpr { sum_data_outer_products_.Init(matrix_dimension, matrix_dimension); data_outer_products_far_field_expansion_.Init(matrix_dimension); + weight_diagram_far_field_expansion_.Init(matrix_dimension); sum_target_weighted_data_.Init(matrix_dimension); target_weighted_data_far_field_expansion_.Init(matrix_dimension); @@ -125,9 +134,11 @@ class DenseLpr { target_weighted_data_far_field_expansion_[j].Init(dimension); data_outer_products_far_field_expansion_[j].Init(matrix_dimension); - + weight_diagram_far_field_expansion_[j].Init(matrix_dimension); + for(index_t i = 0; i < matrix_dimension; i++) { data_outer_products_far_field_expansion_[j][i].Init(dimension); + weight_diagram_far_field_expansion_[j][i].Init(dimension); } } @@ -650,6 +661,10 @@ class DenseLpr { rset_variance_ *= 1.0 / (rset_.n_cols() - 2.0 * rset_first_degree_of_freedom_ + rset_second_degree_of_freedom_); + + printf("First degree of freedom: %g\n", rset_first_degree_of_freedom_); + printf("Second degree of freedom: %g\n", rset_second_degree_of_freedom_); + printf("Reference set variance: %g\n", rset_variance_); } void ComputeConfidenceBands_(const Matrix &queries, @@ -895,9 +910,11 @@ class DenseLpr { stream = fopen(fname, "w+"); } for(index_t r = 0; r < rset_.n_cols(); r++) { - fprintf(stream, "%g %g %g %g\n", rset_confidence_bands_[r].lo, + fprintf(stream, "%g %g %g %g %g %g\n", rset_confidence_bands_[r].lo, rset_regression_estimates_[r], rset_confidence_bands_[r].hi, - leave_one_out_rset_regression_estimates_[r]); + leave_one_out_rset_regression_estimates_[r], + rset_magnitude_weight_diagrams_[r], + rset_influence_values_[r]); } if(stream != stdout) { diff --git a/fastlib2/contrib/dongryel/regression/dense_lpr_impl.h b/fastlib2/contrib/dongryel/regression/dense_lpr_impl.h index 642011d3e4..774f1f51d0 100644 --- a/fastlib2/contrib/dongryel/regression/dense_lpr_impl.h +++ b/fastlib2/contrib/dongryel/regression/dense_lpr_impl.h @@ -110,10 +110,24 @@ void DenseLpr::BasicComputeDualTree_ if(leave_one_out_query_regression_estimates != NULL) { for(index_t i = 0; i < tmp_q_results.length(); i++) { tmp_q_results[old_from_new_queries[i]] = - (*leave_one_out_query_regression_estimates)[i]; + (*leave_one_out_query_regression_estimates)[i]; } leave_one_out_query_regression_estimates->CopyValues(tmp_q_results); } + if(query_magnitude_weight_diagrams != NULL) { + for(index_t i = 0; i < tmp_q_results.length(); i++) { + tmp_q_results[old_from_new_queries[i]] = + (*query_magnitude_weight_diagrams)[i]; + } + query_magnitude_weight_diagrams->CopyValues(tmp_q_results); + } + if(query_influence_values != NULL) { + for(index_t i = 0; i < tmp_q_results.length(); i++) { + tmp_q_results[old_from_new_queries[i]] = + (*query_influence_values)[i]; + } + query_influence_values->CopyValues(tmp_q_results); + } } template @@ -315,6 +329,18 @@ InitializeReferenceStatistics_(ReferenceTree *rnode) { rnode->stat().data_outer_products_far_field_expansion_[j][i]. Add(reference_point_expansion[j] * reference_point_expansion[i] / norm_constant, kernels_[r].bandwidth_sq(), r_col); + + // Note that for computing the magnitude of the weight + // diagram vector, we need to compute the far-field + // expansion of B^T W(q)^2 B, so the thing is a bit tricky + // in the variable-bandwidth case. Note that we divide by + // the squared normalization constant due to the squaring + // the weight matrix W(q). + rnode->stat().weight_diagram_far_field_expansion_[j][i]. + Add(reference_point_expansion[j] * reference_point_expansion[i] / + (norm_constant * norm_constant), kernels_[r].bandwidth_sq(), + r_col); + rnode->stat().sum_data_outer_products_.set (i, j, rnode->stat().sum_data_outer_products_.get(i, j) + reference_point_expansion[j] * @@ -392,6 +418,12 @@ InitializeReferenceStatistics_(ReferenceTree *rnode) { rnode->stat().data_outer_products_far_field_expansion_[j][i]. Add(rnode->right()->stat(). data_outer_products_far_field_expansion_[j][i]); + + rnode->stat().weight_diagram_far_field_expansion_[j][i]. + Add(rnode->left()->stat().weight_diagram_far_field_expansion_[j][i]); + rnode->stat().weight_diagram_far_field_expansion_[j][i]. + Add(rnode->right()->stat(). + weight_diagram_far_field_expansion_[j][i]); } // end of iterating over each row. } // end of iterating over each column. @@ -716,10 +748,10 @@ void DenseLpr::DualtreeLprCanonical_ qnode->stat().postponed_weight_diagram_numerator_l_.set (j, i, qnode->stat().postponed_weight_diagram_numerator_l_.get(j, i) + rnode->stat(). - data_outer_products_far_field_expansion_[j][i]. + weight_diagram_far_field_expansion_[j][i]. ComputeMinKernelSum(qnode->bound())); qnode->stat().postponed_moment_weight_diagram_numerator_e_[j][i]. - Add(rnode->stat().data_outer_products_far_field_expansion_[j][i]); + Add(rnode->stat().weight_diagram_far_field_expansion_[j][i]); } } @@ -1019,9 +1051,33 @@ FinalizeQueryTree_(QueryTree *qnode, const Matrix &qset, (dimension_, lpr_order_, query_point, query_point_expansion.ptr()); (*query_regression_estimates)[q] = la::Dot(query_point_expansion, least_squares_solution); + + // Now we compute the magnitude of the weight diagram for each + // query point. This step depends on the pseudoinverse computed + // above. + Vector pseudo_inverse_times_query_expansion; + Vector intermediate_product; + la::MulInit(pseudoinverse_denominator, query_point_expansion, + &pseudo_inverse_times_query_expansion); + la::MulInit(weight_diagram_numerator_e[q], + pseudo_inverse_times_query_expansion, &intermediate_product); + (*query_magnitude_weight_diagrams)[q] = + sqrt(la::Dot(pseudo_inverse_times_query_expansion, + intermediate_product)); + + // Compute the influence value at each point (if it belongs to + // the reference set), i.e. (r(q))^T (B^T W(q) B)^-1 B^T W(q) + // e_i = (r(q))^T (B^T W(q) B)-1 r(q). + if(query_influence_values != NULL) { + (*query_influence_values)[q] = + la::Dot(query_point_expansion, pseudo_inverse_times_query_expansion); + } // Compute the leave-one-out regression estimate for - // cross-validation stage. + // cross-validation stage. Note that this needs to be after + // computing the magnitude of the weight diagram and the + // influence values, since these two depend on the + // non-leave-one-out quantities. if(leave_one_out_query_regression_estimates != NULL) { // Subtract the contribution of the point itself from the @@ -1052,27 +1108,8 @@ FinalizeQueryTree_(QueryTree *qnode, const Matrix &qset, // query point by the beta_q. (*leave_one_out_query_regression_estimates)[q] = la::Dot(least_squares_solution, query_point_expansion); - } - - // Now we compute the magnitude of the weight diagram for each - // query point. - Vector pseudo_inverse_times_query_expansion; - Vector intermediate_product; - la::MulInit(pseudoinverse_denominator, query_point_expansion, - &pseudo_inverse_times_query_expansion); - la::MulInit(weight_diagram_numerator_e[q], - pseudo_inverse_times_query_expansion, &intermediate_product); - (*query_magnitude_weight_diagrams)[q] = - sqrt(la::Dot(pseudo_inverse_times_query_expansion, - intermediate_product)); + } // end of computing loo regression estimate... - // Compute the influence value at each point (if it belongs to - // the reference set), i.e. (r(q))^T (B^T W(q) B)^-1 B^T W(q) - // e_i = (r(q))^T (B^T W(q) B)-1 r(q). - if(query_influence_values != NULL) { - (*query_influence_values)[q] = - la::Dot(query_point_expansion, pseudo_inverse_times_query_expansion); - } } // end of iterating over each query point. } else { diff --git a/fastlib2/contrib/dongryel/regression/krylov_lpr.h b/fastlib2/contrib/dongryel/regression/krylov_lpr.h index a451305cf6..936a6453d3 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr.h +++ b/fastlib2/contrib/dongryel/regression/krylov_lpr.h @@ -354,7 +354,8 @@ class KrylovLpr { */ void ComputeWeightedVectorSum_ (QueryTree *qroot, const Matrix &qset, const Vector &weights, - const ArrayList *query_in_cg_loop, index_t column_index, + const ArrayList *query_in_cg_loop, + const bool confidence_band_computaton_phase, index_t column_index, Matrix &right_hand_sides_l, Matrix &right_hand_sides_e, Vector &right_hand_sides_used_error, Vector &right_hand_sides_n_pruned, Matrix *leave_one_out_right_hand_sides_e) { @@ -445,7 +446,7 @@ class KrylovLpr { // query point. printf("Starting Phase 1...\n"); ComputeWeightedVectorSum_ - (qroot, qset, rset_target_divided_by_norm_consts_, NULL, 0, + (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); @@ -608,6 +609,11 @@ class KrylovLpr { } } + void LinearOperatorConfidenceBand + (QueryTree *qroot, const Matrix &qset, + const Matrix &query_expansion_solution_vectors, + Matrix &linear_transformed_query_expansion_solution_vectors); + public: ////////// Constructor/Destructor ////////// diff --git a/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h b/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h index 20a17a19b1..3383b12e6d 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h +++ b/fastlib2/contrib/dongryel/regression/krylov_lpr_setup_impl.h @@ -635,7 +635,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)[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 90165f6a4f..5722219ba0 100644 --- a/fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h +++ b/fastlib2/contrib/dongryel/regression/krylov_lpr_solver_impl.h @@ -36,7 +36,7 @@ void KrylovLpr::LinearOperator // Compute the current column linear operator. ComputeWeightedVectorSum_ - (qroot, qset, rset_inv_norm_consts_, &query_in_cg_loop, d, + (qroot, qset, rset_inv_norm_consts_, &query_in_cg_loop, false, d, vector_l, vector_e, vector_used_error, vector_n_pruned, leave_one_out_vector_e); diff --git a/fastlib2/contrib/dongryel/regression/naive_lpr.h b/fastlib2/contrib/dongryel/regression/naive_lpr.h index 73f336ebc0..ab93d81130 100644 --- a/fastlib2/contrib/dongryel/regression/naive_lpr.h +++ b/fastlib2/contrib/dongryel/regression/naive_lpr.h @@ -337,6 +337,8 @@ class NaiveLpr { (rset_.n_cols() - 2.0 * rset_first_degree_of_freedom_ + rset_second_degree_of_freedom_); + printf("First degree of freedom: %g\n", rset_first_degree_of_freedom_); + printf("Second degree of freedom: %g\n", rset_second_degree_of_freedom_); printf("Reference set variance: %g\n", rset_variance_); } @@ -522,9 +524,11 @@ class NaiveLpr { stream = fopen(fname, "w+"); } for(index_t r = 0; r < rset_.n_cols(); r++) { - fprintf(stream, "%g %g %g %g\n", rset_confidence_bands_[r].lo, + fprintf(stream, "%g %g %g %g %g %g\n", rset_confidence_bands_[r].lo, rset_regression_estimates_[r], rset_confidence_bands_[r].hi, - leave_one_out_rset_regression_estimates_[r]); + leave_one_out_rset_regression_estimates_[r], + rset_magnitude_weight_diagrams_[r], + rset_influence_values_[r]); } if(stream != stdout) {