Sketches 1000 times slower in ubuntu than windows

[shoogen]

or use a github mirror https://github.com/tjtg/eigen/compare/3.0.3...3.0.4
a good time to do a "git bisect"

[wmayer]

OK, here we go. The different behaviour is caused by the change inside

Code: Select all

FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)

See https://bitbucket.org/eigen/eigen/diff/ ... at=default

https://bitbucket.org/eigen/eigen/histo ... at=default commit e5fa60b with log message:

stop fill pivoting LU only if the pivot is exactly 0

[wmayer]

Since Eigen3 is a pure template library we can implement our own method FullPivLU::compute() for a special type. So, in this case we can put the following code block which has the old behaviour to GCS.cpp

Code: Select all

namespace Eigen {

template<>
FullPivLU< Matrix<double,-1,-1,0,-1,-1> >& FullPivLU< Matrix<double,-1,-1,0,-1,-1> >::compute(const Matrix<double,-1,-1,0,-1,-1>& matrix)
{
  m_isInitialized = true;
  m_lu = matrix;

  const Index size = matrix.diagonalSize();
  const Index rows = matrix.rows();
  const Index cols = matrix.cols();

  // will store the transpositions, before we accumulate them at the end.
  // can't accumulate on-the-fly because that will be done in reverse order for the rows.
  m_rowsTranspositions.resize(matrix.rows());
  m_colsTranspositions.resize(matrix.cols());
  Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i

  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
  m_maxpivot = RealScalar(0);
  RealScalar cutoff(0);

  for(Index k = 0; k < size; ++k)
  {
    // First, we need to find the pivot.

    // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
    Index row_of_biggest_in_corner, col_of_biggest_in_corner;
    RealScalar biggest_in_corner;
    biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
                        .cwiseAbs()
                        .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
    row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
    col_of_biggest_in_corner += k; // need to add k to them.

    // when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix
    if(k == 0) cutoff = biggest_in_corner * NumTraits<Scalar>::epsilon();

    // if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values.
    // Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in
    // their pseudo-code, results in numerical instability! The cutoff here has been validated
    // by running the unit test 'lu' with many repetitions.
    if(biggest_in_corner < cutoff)
    {
      // before exiting, make sure to initialize the still uninitialized transpositions
      // in a sane state without destroying what we already have.
      m_nonzero_pivots = k;
      for(Index i = k; i < size; ++i)
      {
        m_rowsTranspositions.coeffRef(i) = i;
        m_colsTranspositions.coeffRef(i) = i;
      }
      break;
    }

    if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;

    // Now that we've found the pivot, we need to apply the row/col swaps to
    // bring it to the location (k,k).

    m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
    m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
    if(k != row_of_biggest_in_corner) {
      m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
      ++number_of_transpositions;
    }
    if(k != col_of_biggest_in_corner) {
      m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
      ++number_of_transpositions;
    }

    // Now that the pivot is at the right location, we update the remaining
    // bottom-right corner by Gaussian elimination.

    if(k<rows-1)
      m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
    if(k<size-1)
      m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).tail(rows-k-1) * m_lu.row(k).tail(cols-k-1);
  }

  // the main loop is over, we still have to accumulate the transpositions to find the
  // permutations P and Q

  m_p.setIdentity(rows);
  for(Index k = size-1; k >= 0; --k)
    m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));

  m_q.setIdentity(cols);
  for(Index k = 0; k < size; ++k)
    m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));

  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
  return *this;
}

}

[logari81]

or just replace FullPivLU with a faster version. The only reason for using FullPivLU was that it was supposed to be the most robust and it was extensively tested with many examples (before it was broken).

[triplus]

Unpack them on your disk and use a diff tool to see what has changed.

Just a side note there is nice GUI tool named Meld that works quite good on Linux and it's usually available in official repositories of popular Linux distributions:

http://meldmerge.org/

OK, here we go. The different behaviour is caused by the change inside...

Fantastic then this is probably regression introduced in Eigen 3.0.4? Should we report this to Eigen devs now when we have such detailed information?

Since Eigen3 is a pure template library we can implement our own method FullPivLU::compute() for a special type. So, in this case we can put the following code block which has the old behaviour to GCS.cpp

That would be a quick fix yes and i think it makes most sense if it can be done like that. In the future probably ideal solution would be for Eigen devs to fix the regression and to not having issues anymore by using default Eigen code?

or just replace FullPivLU with a faster version. The only reason for using FullPivLU was that it was supposed to be the most robust and it was extensively tested with many examples (before it was broken).

I am always interested in "fast" but if i understand you correctly current method is used because it should be extensively tested and could changing this to some other method easily introduce new bugs?

Anyway now when you the devs have such detailed info you will know what to do with it... i am just waiting it lands on daily PPA to test it!

[yoshimitsuspeed]

Wow, great to see how much attention this is getting. A lot of it's over my head lol.
Sounds like there is some progress being made though, that's awesome.

[wmayer]

or just replace FullPivLU with a faster version. The only reason for using FullPivLU was that it was supposed to be the most robust and it was extensively tested with many examples

And what do you suggest to use instead? As you said it's the most robust algorithm wouldn't it be best to keep it then? And just to point it our again: the problem of FullPivLU itself is not speed but accuracy which causes solve_DL() to run "forever".

Should we report this to Eigen devs now when we have such detailed information?

Well, that's the question. According to the comment in the old code

// if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values.
// Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in
// their pseudo-code, results in numerical instability! The cutoff here has been validated
// by running the unit test 'lu' with many repetitions.

the algorithm was numerical stable and passed their unit tests. For any reason they changed their mind

stop fill pivoting LU only if the pivot is exactly 0

and followed the pseudo-code they had from a paper. So, maybe we could ask them for their motivation to do this change. But even if they undo this change it doesn't really help us because we need a solution now that works with the Eigen3 packages we have in the repositories of the current distributions.

[triplus]

I see they changed that behaviour intentionally and i am guessing FreeCAD can't just conform to that changes easily? Then alternative solution as logari81 suggested might be the best option because that would fix everything in a way it would be compatible with all Eigen3 packages out there as you said.The only downside is it could bring in additional bugs?

[shoogen]

I see they changed that behaviour intentionally and i am guessing FreeCAD can't just conform to that changes easily?

It might be the case that our matrix is just ill-conditioned. But nevertheless it think they should have been more verbose on their change. i only took one course on numerical mathematics. Having skipped through chapter 3.4.8 3.4.9 of the 3rd edition of "Matrix Computation" by G. Golub and C. van Loan it think the maintainers of Eigen made a mistake. Because though absent in the pseudo code the problem of rounding errors and rank determination is mentioned in the text (esp. in 3.4.9)
I think werner is right. But we should file a bug for eigen, as well.

[wmayer]

Until we have a better solution I checked in my suggestion.