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

This commit is contained in:
Dongryeol Lee
2008-03-02 04:39:00 +00:00
parent 28e2e2b364
commit f6dba35e8b
6 changed files with 97 additions and 33 deletions
@@ -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) {
@@ -110,10 +110,24 @@ void DenseLpr<TKernel, TPruneRule>::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<typename TKernel, typename TPruneRule>
@@ -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<TKernel, TPruneRule>::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 {
@@ -354,7 +354,8 @@ class KrylovLpr {
*/
void ComputeWeightedVectorSum_
(QueryTree *qroot, const Matrix &qset, const Vector &weights,
const ArrayList<bool> *query_in_cg_loop, index_t column_index,
const ArrayList<bool> *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 //////////
@@ -635,7 +635,7 @@ void KrylovLpr<TKernel, TPruneRule>::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;
}
@@ -36,7 +36,7 @@ void KrylovLpr<TKernel, TPruneRule>::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);
@@ -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) {