Last generation of the Krylov subspace based method has been caught - the bug was due to reusing the same dual-tree computation function across three different phases (some phases were required to use squared kernel values, rather than the regular values)

2008-03-03 20:52:57 +00:00
parent 221c1ce793
commit bd9dfdcf9e
5 changed files with 70 additions and 42 deletions
@@ -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<bool> *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_));
    }
  }

@@ -78,7 +78,7 @@ void KrylovLpr<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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<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 != NULL && !((*query_in_cg_loop)[q])) {
 	continue;
      }
      
@@ -30,6 +30,11 @@ void KrylovLpr<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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++) {
@@ -5,9 +5,16 @@
 #endif

 template<typename TKernel, typename TPruneRule>
-void KrylovLpr<TKernel, TPruneRule>::TestRightHandSideComputation_
-(const Matrix &qset, const Matrix &approximated) {
-  
+void KrylovLpr<TKernel, TPruneRule>::TestDualtreeComputation_
+(const Matrix &qset, const ArrayList<bool> *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<TKernel, TPruneRule>::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<TKernel, TPruneRule>::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);
@@ -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");
  }