Hello.

I really stuck with the following problem (least squares).

First I create 2 big matrices:

Code: Select all

int const N = 10000;
Eigen::Matrix<int, 3, Eigen::Dynamic> x(3, N);
Eigen::Matrix<float, 2, Eigen::Dynamic> f(2, N);

Then I fill them by data:

Code: Select all

for(int n(0); n < N; ++n)
{
   int const a( rand() ), b( rand() ), c( rand() );
   float const d(a * 1.1f + b * 2.1f + c * 3.7f + (rand() % 10) - 5);
   float const e(a * 1.5f + b * 2.5f + c * 3.1f + (rand() % 8) - 4);

   x.col(n) << a, b, c;
   f.col(n) << d, e;
}

Then I trying to calculate least squares:

Code: Select all

Eigen::Matrix<double, 3, 3> normalMatrix;
normalMatrix.triangularView<Eigen::Upper>() = x.cast<double>() * x.cast<double>().transpose(); //OK

Eigen::Matrix<double, 3, 2> dataMatrix;
dataMatrix.noalias() = x.cast<double>() * f.cast<double>().transpose(); //STACK OVERFLOW HERE

Eigen::Matrix<float, 2, 3> result;
result = Eigen::LDLT< Eigen::Matrix<double, 3, 3> >( normalMatrix.selfadjointView<Eigen::Upper>() ).solve(dataMatrix).transpose().cast<float>();

I cannot figure out why I can produce 3×3 matrix by doing x.cast<double>() * x.cast<double>().transpose(), but cannot produce 3×2 matrix by doing x.cast<double>() * f.cast<double>().transpose().

I am sure that I can substitute for-loops doing multiplication instead of calling eigen multiplication, and program will work.

P.S. The type conversions are necessary to me. I need to find matrix that maps 3-component integer vectors into 2-components float vectors. The values of matrix x * x.transpose() will not fit into integer, so I want to use double to do summation.

If I use N = 1000 insted of N = 10000 then I get correct result:

Code: Select all

std::cout << result;

Gives me

Code: Select all

1.09998     2.1 3.69999
1.49999     2.5 3.09998

But I need to use millions of datapoints.

sorry but I cannot reproduce.

Can you give us information about your system and compiler.

here is a complete example:

Code: Select all

#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;

int main()
{
  int const N = 100000;
  Eigen::Matrix<int, 3, Eigen::Dynamic> x(3, N);
  Eigen::Matrix<float, 2, Eigen::Dynamic> f(2, N);
 
  for(int n(0); n < N; ++n)
  {
    int const a( rand() ), b( rand() ), c( rand() );
    float const d(a * 1.1f + b * 2.1f + c * 3.7f + (rand() % 10) - 5);
    float const e(a * 1.5f + b * 2.5f + c * 3.1f + (rand() % 8) - 4);

    x.col(n) << a, b, c;
    f.col(n) << d, e;
  }
 
  Eigen::Matrix<double, 3, 3> normalMatrix;
  normalMatrix.triangularView<Eigen::Upper>() = x.cast<double>() * x.cast<double>().transpose(); //OK

  Eigen::Matrix<double, 3, 2> dataMatrix;
  dataMatrix.noalias() = x.cast<double>() * f.cast<double>().transpose(); //STACK OVERFLOW HERE

  Eigen::Matrix<float, 2, 3> result;
  result = Eigen::LDLT< Eigen::Matrix<double, 3, 3> >( normalMatrix.selfadjointView<Eigen::Upper>() ).solve(dataMatrix).transpose().cast<float>();

  std::cout << result << "\n";
 
  return 0;
}



Interestingly that your code works fine, while my code crashes. The only difference of my code from yours is presence of templated function (I do not know why this matters):

Code: Select all

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Cholesky>

template<class CALC, int N, int M, class TX, class TF, class TR, int P> static inline void solveLDLT(Eigen::Matrix<TX, N, P> const &x, Eigen::Matrix<TF, M, P> const &f, Eigen::Matrix<TR, M, N> &result)
{
   typedef Eigen::Matrix<CALC, N, N> NormalMatrix;
   typedef Eigen::Matrix<CALC, N, M> DataMatrix;

   NormalMatrix normalMatrix;
   normalMatrix.triangularView<Eigen::Upper>() = x.cast<CALC>() * x.cast<CALC>().transpose();
   DataMatrix dataMatrix;
   dataMatrix.noalias() = x.cast<CALC>() * f.cast<CALC>().transpose();
   result = Eigen::LDLT<NormalMatrix>( normalMatrix.selfadjointView<Eigen::Upper>() ).solve(dataMatrix).transpose().cast<TR>();
}

int main(int argc, char **argv)
{
   int const N = 10000;
   Eigen::Matrix<int, 3, Eigen::Dynamic> x(3, N);
   Eigen::Matrix<float, 2, Eigen::Dynamic> f(2, N);
   
   for(int n(0); n < N; ++n)
   {
      int const a( rand() ), b( rand() ), c( rand() );
      float const d(a * 1.1f + b * 2.1f + c * 3.7f + (rand() % 10) - 5);
      float const e(a * 1.5f + b * 2.5f + c * 3.1f + (rand() %  8) - 4);

      x.col(n) << a, b, c;
      f.col(n) << d, e;
   }

   Eigen::Matrix<double, 2, 3> result;

   solveLDLT<double, 3, 2, int, float, double, N>(x, f, result);
   std::cout << result << std::endl << std::endl;

   return 0;
}

I use Visual Studio 2010 Express (crash on both Debug and Release default settings), Windows Vista (32-bit).

Sorry! Now I see the reason for the error!

I should have used Eigen::Dynamic instead of "N" as the last template argument!