I am trying to calculate x.transpose()*x where x is a vector in two different ways and I don't get the same answer. My code is below

Code: Select all

#include <cstdlib>
#include <cmath>
#include <limits>
#include <iostream>
#include <Eigen/Core>

int test_eigen_product(void)
{
   typedef double RealType;
   typedef Eigen::Matrix<RealType, Eigen::Dynamic, 1> RealVectorType;
   typedef RealVectorType::Index IndexType;
   
   const IndexType N=1000000;
   RealVectorType x=RealVectorType::Random(N);
   RealType val=0;
   
   val=x.transpose()*x;
   
   RealType valTest=0;
   for(IndexType i=0;i<N;++i)
   {
      valTest+=x(i)*x(i);
   }
   
   RealType eps=std::numeric_limits<RealType>::epsilon();
   
   if(fabs(val-valTest) > eps )
   {
      std::cout<<"fabs(val-valTest)= "<<fabs(val-valTest)<<std::endl;
      return EXIT_FAILURE;
   }
   
   return EXIT_SUCCESS;
}

int main(void)
{
   int ret=EXIT_SUCCESS;
   ret+=test_eigen_product();
   return ret;
}



I get fabs(val-valTest)= 5.82077e-09. Is the difference coming from the addition valTest+=x(i)*x(i) ? Many thanks for your help.

Compiler is gcc 4.8.1 with flags -std=c++11 -pedantic -Wall -Wextra -Wfatal-errors -g

Sree

This seems to be just the normal roundoff error. You are adding a million numbers, randomly distributed between 0 and 1, so the sum is approximately 500,000. Thus an error of 5e-9 translates to a relative error of 5e-14 or about 200 times epsilon, which seems reasonable for a million operations. If you compare the loop that you wrote by hand with the same loop backwards (as below), you get a similar difference in the result.

Code: Select all

for(IndexType i=N-1; i>=0; --i)
{
   valTest+=x(i)*x(i);
}


Makes sense. I ran the loop you suggested to see the result. I got fabs(val-valTest2)= 9.02219e-09 which is comparable to earlier value. Thank you.
Sree