Changeset 4053

Timestamp:

Mar 26, 2021, 4:29:42 AM (22 months ago)

Author:

Peter

Message:

add some test for correlation(Matrix); speed up, particuarly when #rows is much greater than #columns. Avoid copying and modifying the whole matrix, instead calculate simple sums and squared sums using BLAS functionality as much as possible.

Location:

trunk

Files:

• test/statistics.cc (modified) (2 diffs)
• yat/statistics/utility.cc (modified) (1 diff)
• yat/statistics/utility.h (modified) (1 diff)

trunk/test/statistics.cc

@@ -150 +150 @@
 void test_correlation(test::Suite& suite)
 {
-  utility::Matrix x(10,100);
+  utility::Matrix x(100,10);
+  for (size_t i=0; i<x.rows(); ++i)
+    for (size_t j=0; j<x.columns(); ++j)
+      x(i, j) = sin(5*i + j);
   utility::Matrix C = statistics::correlation(x);
   if (!suite.equal(C(0,0), 1.0)) {
@@ -156 +159 @@
     suite.err() << "correlation element(0, 0) not unity\n";
   }
+
+  for (size_t i=0; i<C.rows(); ++i)
+    for (size_t j=0; j<C.columns(); ++j) {
+      statistics::AveragerPair ap;
+      add(ap, x.begin_column(i), x.end_column(i), x.begin_column(j));
+      if (!suite.equal(C(i, j), ap.correlation(), 10)) {
+        suite.add(false);
+        suite.err() << "error: C(" << i << ", " << j << ") is " << C(i,j)
+                    << "; expected " << ap.correlation() << "\n";
+      }
+    }
+
 }
 

trunk/yat/statistics/utility.cc

@@ -50 +50 @@
 
 
-  utility::Matrix correlation(utility::Matrix x)
+  utility::Matrix correlation(const utility::Matrix& X)
   {
-    using utility::Matrix;
-    size_t n = x.columns();
-    // centralise
-    normalizer::ColumnNormalizer<normalizer::Centralizer<> > normalizer;
-    normalizer(x, x);
-
-    Matrix cov(n, n);
-    for (size_t i=0; i<n; ++i) {
-      for (size_t j=i; j<n; ++j) {
-        cov(i, j) = x.column_const_view(i) * x.column_const_view(j);
-        cov(j, i) = cov(i, j);
-      }
+    size_t cols = X.columns();
+    size_t rows = X.rows();
+    utility::Vector m(cols);
+    utility::Vector x2(cols);
+    utility::Vector stddev(cols);
+    for (size_t i=0; i<cols; ++i) {
+      utility::VectorConstView vec = X.column_const_view(i);
+      m(i) = sum(vec);
+      x2(i) = vec * vec;
+      // scaled standard deviation
+      stddev(i) = std::sqrt(x2(i) - m(i) * m(i) / rows);
     }
 
-    utility::Vector stddev(n);
-    for (size_t i=0; i<n; ++i)
-      stddev(i) = std::sqrt(cov(i, i));
+    utility::Matrix corr(cols, cols, 1.0);
+    for (size_t i=0; i<cols; ++i)
+      for (size_t j=i+1; j<cols; ++j) {
+        corr(i,j) =
+          X.column_const_view(i) * X.column_const_view(j) - m(i)*m(j)/rows;
+        corr(i, j) /= stddev(i) * stddev(j);
+        // symmetry
+        corr(j, i) = corr(i,j);
+      }
 
-    Matrix corr(cov);
-    for (size_t i=0; i<n; ++i) {
-      corr(i, i) = 1.0;
-      for (size_t j=0; j<i; ++j) {
-        corr(i, j) = cov(i, j) / (stddev(i) * stddev(j));
-        corr(j, i) = corr(i, j);
-      }
-    }
     return corr;
   }

trunk/yat/statistics/utility.h

@@ -145 +145 @@
      \since New in yat 0.16
    */
-  utility::Matrix correlation(utility::Matrix x);
+  utility::Matrix correlation(const utility::Matrix& x);