I have two programs. Both perform the same calculation, but have different results.
Can someone help interpret why the results are different?

Program tmp.cpp gets the result of -23.8359, and the core calculation is:

Code: Select all

  Eigen::MatrixXf ret = mat.topRows(31795).array().square().matrix().colwise().sum() -
      mat.bottomRows(1).array().square().matrix().colwise().sum();

Program tmp2.cpp gets the result of 11.3672 , and the core calculation is:

Code: Select all

  Eigen::MatrixXf m = mat.topRows(31795).array().square();
  Eigen::VectorXf v = m.colwise().sum();
      m = mat.bottomRows(1).array().square();
    v -= m.colwise().sum();

I have put all source codes below, and the data file "tmp", a column data matrix (31796 by 1), to this URL http://s000.tinyupload.com/index.php?file_id=04938324134919232586.
My environment is gcc (GCC) 7.1.0, Eigen 3.3.4.

tmp.cpp

Code: Select all

#include <iostream>
#include <Eigen/Dense>
using namespace std;
int main()
{
  Eigen::MatrixXf mat(31796,1);
  FILE* fp = fopen("tmp", "r");
  float a;
  for (int i = 0; i < 31796; ++i) {
     fscanf(fp, "%f", &a);
     mat(i, 0) = a;
  }

  Eigen::MatrixXf ret = mat.topRows(31795).array().square().matrix().colwise().sum() -
      mat.bottomRows(1).array().square().matrix().colwise().sum();

  std::cout << ret << std::endl;
}


tmp2.cpp

Code: Select all

#include <iostream>
#include <Eigen/Dense>
using namespace std;
int main()
{
  Eigen::MatrixXf mat(31796,1);
  FILE* fp = fopen("tmp", "r");
  float a;
  for (int i = 0; i < 31796; ++i) {
     fscanf(fp, "%f", &a);
     mat(i, 0) = a;
  }

  Eigen::MatrixXf m = mat.topRows(31795).array().square();
  Eigen::VectorXf v = m.colwise().sum();
      m = mat.bottomRows(1).array().square();
    v -= m.colwise().sum();

  std::cout << v << std::endl;
}


Since both codes do not performs the arithmetic operations in the exact same order, because of rounding errors it is expected to get slightly different results. In your case (i.e., for your precise computation AND numerical values), this is worst because you hit catastrophic cancellation issues (when doing the '-' operation), so you need to increase the numerical precision to prevent them, for instance using MatrixXd and VectorXd, and then you get 13.2306 in both cases. You can also check that the single precision codes are correct by replacing the numerical values of mat by:

mat.random();

and then, you'll see that the core of the problem comes from your particular numerical values.

BTW, better use VectorX? instead of MatrixX? whenever you know you are dealing with vectors.

Hi ggael,

I am not sure why the codes "do not performs the arithmetic operations in the exact same order".

To make sure codes performs the arithmetic operations in the exact same order and catastrophic cancellation will not happen during subtraction, I revised both programs.
Now there are only square and sum operations: both will square the first 31795 elements and sum them up.

But the results are still not the same.

BTW: the first 31795 elements include 0, 1, 1.94737 and 2.

PS: I can use VectorX instead of MatrixX in this toy example, but my real application has to use MatrixX.
PS2: I can see

Code: Select all

mat = Eigen::MatrixXf::Random()

works. But I still do not understand why data set will not work.

tmp3.cpp
result = 120549

Code: Select all

#include <iostream>
#include <Eigen/Dense>
using namespace std;
int main()
{
  Eigen::MatrixXf mat(31796,1);
  FILE* fp = fopen("tmp", "r");
  float a;
  for (int i = 0; i < 31796; ++i) {
     fscanf(fp, "%f", &a);
     mat(i, 0) = a;
  }

  Eigen::MatrixXf ret = mat.topRows(31795).array().square().matrix().colwise().sum();

  std::cout << ret << std::endl;
}

tmp4.cpp
result = 120584

Code: Select all

#include <iostream>
#include <Eigen/Dense>
using namespace std;
int main()
{
  Eigen::MatrixXf mat(31796,1);
  FILE* fp = fopen("tmp", "r");
  float a;
  for (int i = 0; i < 31796; ++i) {
     fscanf(fp, "%f", &a);
     mat(i, 0) = a;
  }

  Eigen::MatrixXf m = mat.topRows(31795).array().square();
  Eigen::VectorXf v = m.colwise().sum();

  std::cout << v << std::endl;
}

Did you suggest the internal computation are different for these two statements?

Code: Select all

  Eigen::MatrixXf ret = mat.topRows(31795).array().square().matrix().colwise().sum();


Code: Select all

  Eigen::MatrixXf m = mat.topRows(31795).array().square();
  Eigen::VectorXf v = m.colwise().sum();


To my knowledge, the latter just split one statement into two statements, so they seem to be equivalent.

Do you think Eigen treats them differently? If so, is there a way I can examine its internal calculations?

Thanks.

You are right, in this case there is no reason that different code to be generated. The problem was a too conservative rule to enable vectorization in Block<> expression. This is fixed in devel and 3.3 branch. Nonetheless, regarding your original issue, I would still recommend switching to double.

BTW, it simpler to write:

Code: Select all

MatrixXf ret = mat.topRows(31795).colwise().squaredNorm() -  mat.bottomRows(1).colwise().squaredNorm();