0

This problem is related to Workaround for resizing Eigen::Ref, however, I do not have the restriction of trying to avoid templates (in-fact, I would like to have a solution working with templates)

I'm using the the eigen library (version 3.2.9, but testet with the latest version of eigen as well with the same outcome) for some automatic differentiation (AD) calculations and came across this "bug", which is already known to the developers (see also the following bug reports: Bug report 1, Bug report 2 and Bug report 3). TL;DR: It is not really a bug, but a generic and clean workaround would likely require some extensive work and since it is not supported by Eigen, not pursued (I guess ...). For me, I am only interested in a subset for of expression for which to get it to work (at least for now) for which there might be a acceptable workaround.

The problem is the following, consider this simplified code, where we have a fixed and dynamic matrix sized AD type

#include <Eigen/Core>
#include "unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h"

typedef Eigen::AutoDiffScalar<Eigen::Matrix<double, -1, 1>> T_dynamic;
typedef Eigen::AutoDiffScalar<Eigen::Matrix<double, 40, 1>> T_fixed;

int main() {
  const T_fixed   fixed   = 10.0;
  const T_dynamic dynamic = 100.0;

  const auto result = fixed - dynamic;

  return 0;
}

This works great. When we hit fixed - dynamic, Eigen will call first the overloaded operator- which looks at both left and right hand side and resizes the derivatives accordingly through the make_coherent_impl() method, which is a template specialised version of make_coherent()

template<typename OtherDerType>
inline const AutoDiffScalar<CwiseBinaryOp<internal::scalar_difference_op<Scalar>, const DerType,const typename internal::remove_all<OtherDerType>::type> >
operator-(const AutoDiffScalar<OtherDerType>& other) const
{
  internal::make_coherent(m_derivatives, other.derivatives());
  return AutoDiffScalar<CwiseBinaryOp<internal::scalar_difference_op<Scalar>, const DerType,const typename internal::remove_all<OtherDerType>::type> >(
  m_value - other.value(),
  m_derivatives - other.derivatives());
}

// resize a to match b is a.size()==0, and conversely.
template<typename A, typename B>
void make_coherent(const A& a, const B&b)
{
  make_coherent_impl<A,B>::run(a.const_cast_derived(), b.const_cast_derived());
}

template<typename A_Scalar, int A_Rows, int A_Cols, int A_Options, int A_MaxRows, int A_MaxCols, typename B>
struct make_coherent_impl<Matrix<A_Scalar, A_Rows, A_Cols, A_Options, A_MaxRows, A_MaxCols>, B> {
  typedef Matrix<A_Scalar, A_Rows, A_Cols, A_Options, A_MaxRows, A_MaxCols> A;
  static void run(A& a, B& b) {
    if((A_Rows==Dynamic || A_Cols==Dynamic) && (a.size()==0))
    {
      a.resize(b.size());
      a.setZero();
    }
  }
};

This works, because both types A and B are of Eigen::Matrix<...> type (there are other permutations of the make_coherent_impl() version which is omitted here).

However, consider the following slightly modified example:

#include <Eigen/Core>
#include "unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h"

typedef Eigen::AutoDiffScalar<Eigen::Matrix<double, -1, 1>> T_dynamic;
typedef Eigen::AutoDiffScalar<Eigen::Matrix<double, 40, 1>> T_fixed;

int main() {

  const double    scalar  = 1.0;
  const T_fixed   fixed   = 10.0;
  const T_dynamic dynamic = 100.0;

  const auto result = dynamic * scalar - fixed * scalar;

  return 0;
}

now, when we hit dynamic * scalar - fixed * scalar and go through the the make_coherent() method, we call

template<typename A, typename B>
struct make_coherent_impl {
  static void run(A&, B&) {}
};

instead, as A and B are no longer of type Eigen::matrix<...> but rather a CwiseUnaryOp<Operation, DerivativeType>.

This type is again a matrix in the end and I want to resize it, but when doing so, eigen detects that it is not a matrix type and calls a no-op resize function (in the DenseBase.h), as it only allows to resize matrices and arrays (see also the function description, here the CwiseUnaryOp is an expression)

/** Only plain matrices/arrays, not expressions, may be resized; therefore the only useful resize methods are
* Matrix::resize() and Array::resize(). The present method only asserts that the new size equals the old size, and does
* nothing else.
*/
void resize(Index newSize)
{
  EIGEN_ONLY_USED_FOR_DEBUG(newSize);
  eigen_assert(newSize == this->size() && "DenseBase::resize() does not actually allow to resize.");
}

The partially specialised template version I am using to call the resize function is as follows

template<class UnaryOpLhs, class DerTypeLhs, class UnaryOpRhs, class DerTypeRhs>
struct make_coherent_impl<Eigen::CwiseUnaryOp<UnaryOpLhs, DerTypeLhs>, Eigen::CwiseUnaryOp<UnaryOpRhs, DerTypeRhs> >
{
  typedef Eigen::CwiseUnaryOp<UnaryOpLhs, DerTypeLhs> A;
  typedef Eigen::CwiseUnaryOp<UnaryOpRhs, DerTypeRhs> B;

  static void run(A& a, B& b) {
    if(a.size()==0 && b.size()!=0)
      a.resize(b.size());
    else if(b.size()==0 && a.size()!=0)
      b.resize(a.size());
  }
};

Is there a way to cast the CwiseUnaryOp to a matrix so we can resize it again or a different route that would achieve the same thing here? I am only showing CwiseUnaryOp here but it should work equally for CwiseBinaryOp

tom
  • 361
  • 3
  • 11
  • Thanks for finding a few duplicates in the bug-tracker :) Do you have use-cases where you actually need to store scalars (without derivatives) inside AD-types? If not, the easiest workaround would be to never store them as AD types (maybe difficult in generic code). The bug indeed has no high priority since it is in an "unsupported" module (proposals for fixes are still welcome). – chtz Nov 07 '19 at 16:53
  • well, the problem is that we don't want to change how we write our code that depends on the AD types (I have suggested ways of how to work around these issues by first resizing all dynamic AD types which works perfectly fine, but that leads to ugly code and in a way is quite restrictive). Without wanting to go into details, the actual code (not this simplification), requires the exact same structure as this example (the substraction of two CwiseUnaryOp of double with Eigen::Matrix) so there is a use case for me ... is it _safe_ to overload the resize method in DenseBase for CwiseUnaryOps? – tom Nov 07 '19 at 20:24

0 Answers0