added correlation matrix computation, including for vd-cma, ref #129

src/cmasolutions.cc

@@ -125,6 +125,33 @@ namespace libcmaes
       }
   }
 
+  dMat CMASolutions::corr() const
+  {
+    dMat corr, dinvcov;
+    if (_cov.size() > 0) // full cov
+      {
+	dinvcov = _cov.diagonal().cwiseSqrt().cwiseInverse();
+	corr = dMat(_cov.rows(),_cov.cols());
+	for (int i=0;i<_cov.cols();i++)
+	  corr.col(i) = _cov.col(i).cwiseProduct(dinvcov);
+	for (int i=0;i<_cov.rows();i++)
+	  corr.row(i) = corr.row(i).cwiseProduct(dinvcov.transpose());
+      }
+    else if (_v.size() > 0) //vd
+      {
+	// we need to compute the full covariance matrix, which is counter productive in large-scale settings
+	dMat cov =_sepcov.asDiagonal()*(dMat::Identity(_sepcov.rows(),_sepcov.rows())+_v*_v.transpose())*(_sepcov.asDiagonal());
+	dinvcov = cov.diagonal().cwiseSqrt().cwiseInverse();
+	corr = dMat(cov.rows(),cov.cols());
+	for (int i=0;i<cov.cols();i++)
+	  corr.col(i) = cov.col(i).cwiseProduct(dinvcov);
+	for (int i=0;i<cov.rows();i++)
+	  corr.row(i) = corr.row(i).cwiseProduct(dinvcov.transpose());
+      }
+    else return dMat::Constant(_sepcov.rows(),1,1.0); // sep
+    return corr;
+  }
+
   void CMASolutions::update_eigenv(const dVec &eigenvalues,
 				   const dMat &eigenvectors)
   {

src/cmasolutions.h

@@ -293,6 +293,12 @@ namespace libcmaes
       return dVec(); // should never reach here.
     }
 
+    /**
+     * \brief returns correlation matrix
+     * @return correlation matrix
+     */
+    dMat corr() const;
+
     /**
      * \brief returns current value of step-size sigma
      * @return current step-size

tests/test-functions.cc

@@ -484,6 +484,8 @@ DEFINE_int32(restarts,9,"maximum number of restarts, applies to IPOP and BIPOP a
 DEFINE_bool(with_gradient,false,"whether to use the function gradient when available in closed form");
 DEFINE_bool(with_num_gradient,false,"whether to use numerical gradient injection");
 DEFINE_bool(with_edm,false,"whether to compute expected distance to minimum when optimization has completed");
+DEFINE_bool(with_stds,false,"whether to compute and print the standard deviation for every parameter");
+DEFINE_bool(with_corr,false,"whether to compute and print the correlation matrix (may not fit in memory in large-scale settings)");
 DEFINE_bool(mt,false,"whether to use parallel evaluation of objective function");
 DEFINE_bool(initial_fvalue,false,"whether to compute initial objective function value at x0");
 DEFINE_int32(elitist,0,"whether to activate elistism, 0: deactivated, 1: reinjects best seen candidate, 2: initial elitism, reinjects x0, 3: on restart scheme, useful when optimizer appears to converge to a value that is higher than the best value reported along the way");
@@ -626,6 +628,10 @@ CMASolutions cmaes_opt()
     }
   if (FLAGS_with_edm)
     LOG(INFO) << "EDM=" << cmasols.edm() << " / EDM/fm=" << cmasols.edm() / cmasols.best_candidate().get_fvalue() << std::endl;
+  if (FLAGS_with_stds)
+    std::cout << "stds=" << cmasols.stds(cmaparams).transpose() << std::endl;
+  if (FLAGS_with_corr)
+    std::cout << "correlation=" << cmasols.corr() << std::endl;
   return cmasols;
 }