Another fix to the naive lpr code

fastlib2/contrib/dongryel/regression/dense_lpr.h

@@ -476,6 +476,8 @@ class DenseLpr {
 
   public:
   
+    ////////// Constructor/Destructor //////////
+  
     /** @brief The constructor which sets pointers to NULL. */
     DenseLpr() {
       qroot_ = NULL;
@@ -493,6 +495,21 @@ class DenseLpr {
       }
     }
 
+    ////////// Getter/Setters //////////
+  
+    /** @brief Get the regression estimates.
+     *
+     *  @param results The uninitialized vector which will be filled
+     *                 with the computed regression estimates.
+     */
+    void get_regression_estimates(Vector *results) { 
+      results->Init(regression_estimates_.length());
+      
+      for(index_t i = 0; i < regression_estimates_.length(); i++) {
+	(*results)[i] = regression_estimates_[i];
+      }
+    }
+
     /////////// User-level Functions //////////
 
     void Compute() {

fastlib2/contrib/dongryel/regression/dense_lpr_main.cc

@@ -1,5 +1,6 @@
 #include "mlpack/kde/dataset_scaler.h"
 #include "dense_lpr.h"
+#include "naive_lpr.h"
 #include "relative_prune_lpr.h"
 
 int main(int argc, char *argv[]) {
@@ -48,12 +49,24 @@ int main(int argc, char *argv[]) {
   DatasetScaler::TranslateDataByMin(queries, references, false);
 
   // Declare local linear krylov object.
-  DenseLpr<GaussianKernel, 0, RelativePruneLpr> local_linear;
-  local_linear.Init(queries, references, reference_targets, 
-		    local_linear_module);
-  local_linear.Compute();
-  local_linear.PrintDebug();
-  
+  Vector fast_lpr_results;
+  DenseLpr<GaussianKernel, 0, RelativePruneLpr> fast_lpr;
+  fast_lpr.Init(queries, references, reference_targets, local_linear_module);
+  fast_lpr.Compute();
+  fast_lpr.PrintDebug();
+  fast_lpr.get_regression_estimates(&fast_lpr_results);
+
+  // Do naive algorithm.
+  Vector naive_lpr_results;
+  NaiveLpr<GaussianKernel, 0> naive_lpr;
+  naive_lpr.Init(queries, references, reference_targets, local_linear_module);
+  naive_lpr.Compute();
+  naive_lpr.PrintDebug();
+  naive_lpr.get_regression_estimates(&naive_lpr_results);
+  printf("Maximum relative error: %g\n", 
+	 MatrixUtil::MaxRelativeDifference(naive_lpr_results, 
+					   fast_lpr_results));
+ 
   // Finalize FastExec and print output results.
   fx_done();
   return 0;

fastlib2/contrib/dongryel/regression/matrix_util.h

@@ -112,6 +112,33 @@ class MatrixUtil {
       }
       return norm_diff;
     }
+
+    static double MaxRelativeDifference(const Vector &true_results,
+					const Vector &approx_results) {
+
+      double max_relative_error = 0;
+
+      for(index_t d = 0; d < true_results.length(); d++) {
+
+	if(isnan(approx_results[d]) || isinf(approx_results[d]) ||
+	   isnan(true_results[d]) || isinf(true_results[d])) {
+	  printf("Warning: Got infinites and NaNs!\n");
+	}
+
+	max_relative_error =
+	  std::max(max_relative_error, 
+		   (fabs(approx_results[d]) - fabs(true_results[d])) /
+		   fabs(true_results[d]));
+
+	printf("%g against %g gives %g\n",
+	       approx_results[d], true_results[d],
+	       (fabs(approx_results[d]) - fabs(true_results[d])) /
+	       fabs(true_results[d]));
+      }
+
+      return max_relative_error;
+    }
+  
 };
 
 #endif

fastlib2/contrib/dongryel/regression/naive_lpr.h

@@ -12,7 +12,7 @@
 #include "multi_index_util.h"
 #include "fastlib/fastlib.h"
 
-template<typename TKernel, int lpr_order = 1>
+template<typename TKernel, int lpr_order>
 class NaiveLpr {
 
   FORBID_ACCIDENTAL_COPIES(NaiveLpr);
@@ -61,7 +61,9 @@ class NaiveLpr {
   int dimension_;
 
  public:
-  
+
+  ////////// Constructor/Destructor //////////
+
   /** @brief The constructor which does nothing.
    */
   NaiveLpr() {}
@@ -70,15 +72,30 @@ class NaiveLpr {
    */
   ~NaiveLpr() {}
 
+  ////////// Getter/Setters //////////
+  
+  /** @brief Get the regression estimates.
+   *
+   *  @param results The uninitialized vector which will be filled
+   *                 with the computed regression estimates.
+   */
+  void get_regression_estimates(Vector *results) { 
+    results->Init(regression_values_.length());
+    
+    for(index_t i = 0; i < regression_values_.length(); i++) {
+      (*results)[i] = regression_values_[i];
+    }
+  }
+
   /** @brief Compute the local polynomial regression values using the
    *         brute-force algorithm.
    */
   void Compute() {
 
     // Temporary variable for storing multivariate expansion of a
-    // reference point.
-    Vector reference_point_expansion;
-    reference_point_expansion.Init(total_num_coeffs_);
+    // point.
+    Vector point_expansion;
+    point_expansion.Init(total_num_coeffs_);
 
     printf("\nStarting naive local polynomial of order %d...\n", lpr_order);
     fx_timer_start(NULL, "naive_local_linear_compute");
@@ -98,7 +115,7 @@ class NaiveLpr {
 
 	// Compute the reference point expansion.
 	MultiIndexUtil::ComputePointMultivariatePolynomial
-	  (dimension_, lpr_order, r_col, reference_point_expansion.ptr());
+	  (dimension_, lpr_order, r_col, point_expansion.ptr());
 	
 	// Compute the pairwise distance and the resulting kernel value.
 	double dsqd = la::DistanceSqEuclidean(qset_.n_rows(), q_col, r_col);
@@ -106,14 +123,13 @@ class NaiveLpr {
 
 	for(index_t i = 0; i < total_num_coeffs_; i++) {
 
-	  numerator_[q][i] += r_target * kernel_value * 
-	    reference_point_expansion[i];
+	  numerator_[q][i] += r_target * kernel_value * point_expansion[i];
 	  
 	  // Here, compute each component of the denominator matrix.
 	  for(index_t j = 0; j < total_num_coeffs_; j++) {
 	    denominator_[q].set(j, i, denominator_[q].get(j, i) +
-				reference_point_expansion[j] * 
-				reference_point_expansion[i] * kernel_value);
+				point_expansion[j] * point_expansion[i] * 
+				kernel_value);
 	  } // End of looping over each (j, i)-th component of the
 	    // denominator matrix.
 	} // End of looping over each i-th component of the numerator
@@ -131,6 +147,10 @@ class NaiveLpr {
       const double *q_col = qset_.GetColumnPtr(q);
       Vector beta_q;
       
+      // Compute the query point expansion.
+      MultiIndexUtil::ComputePointMultivariatePolynomial
+	(dimension_, lpr_order, q_col, point_expansion.ptr());
+
       // Now invert the denominator matrix for each query point and
       // multiply by the numerator vector.
       MatrixUtil::PseudoInverse(denominator_[q], &denominator_inv_q);      
@@ -138,10 +158,7 @@ class NaiveLpr {
 
       // Compute the dot product between the multiindex vector for the
       // query point by the beta_q.
-      regression_values_[q] = beta_q[0];
-      for(index_t i = 1; i <= qset_.n_rows(); i++) {
-	regression_values_[q] += beta_q[i] * q_col[i - 1];
-      }
+      regression_values_[q] = la::Dot(beta_q, point_expansion);
     }
     
     fx_timer_stop(NULL, "naive_local_linear_compute");