Hi,

X is a matrix and y is the column vector corresponding to X's righmost column.

I want to compute the product "yT X", which I call sumXy in the code. There are two ways to compute this with libeigen. One way is just "y.transpose() * x". The other way is to compute each row's product with x, accumulating the result, all done using libeigen.

I find that the two ways don't always give the same exact result. Is this expected? If not, I have a bug. Otherwise, what's one way libeigen could be computing "y.transpose() * x" that results in such a difference? I'm running this on x86-64.

Code: Select all

#include <cassert>
#include <iostream>
#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/Sparse>

using namespace Eigen;

namespace {
double const kCoef = 0.1;
}

bool test(size_t numRows, size_t numCols) {
    MatrixXd x(numRows, numCols);
    for (size_t row = 0; row < numRows; ++row)
        for (size_t col = 0; col < numCols; ++col)
            x(row, col) = (row * 20 + col) * kCoef;
    VectorXd y = x.col(numCols - 1);

    VectorXd sumXy(numCols);
    sumXy.setZero();
    for (size_t row = 0; row < numRows; ++row)
        sumXy += x.row(row) * y[row];
    return sumXy.transpose() == y.transpose() * x;
}

int main() {
    for (size_t numRows = 1; numRows < 20; ++numRows)
        for (size_t numCols = 1; numCols < 20; ++numCols)
            if (!test(numRows, numCols))
                std::cerr << "Failed for (" << numRows << ", " << numCols << ")\n";
    return 0;
}


yes this has to be expected with floating point precision because, for instance, (a+b)+c is different than a+(b+c). In your case, the difference should be in the order of 1e-14. You can check the relative error as follows:
[code]
sumXy2 = y.transpose() * x;
cout << (sumXy2-sumXy).norm() / sumXy2.norm() << "\n";
[code]