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) {