How to keep track of Eigen objects through a solve()?

Question

This question is related to cast from Eigen::CwiseBinaryOp to MatrixXd causes segfault . It will probably have as simple a solution as the former.

In this minimal example, I define Holder, which holds and Eigen matrix, and returns it via its get() member function. Similarly, Decomp is an expression template for the LDLT decomposition of this matrix, and Solve solves for AX=B, yielding X.

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

template <class EigenType> class Holder {
public:
typedef EigenType result_type;

private:
result_type in_;

public:
Holder(const EigenType& in) : in_(in) {}
const result_type& get() const { return in_; }
};

template <class Hold> class Decomp {
public:
typedef typename Eigen::LDLT
    <typename Eigen::MatrixBase<typename Hold::result_type>::PlainObject>
        result_type;

private:
Hold mat_;

public:
Decomp(const Hold& mat) : mat_(mat) {}

result_type get() const { return mat_.get().ldlt(); }
};

template <class Derived, class OtherDerived> class Solve {
public:
typedef typename Eigen::internal::solve_retval
    <typename Derived::result_type, typename OtherDerived::result_type>
        result_type;

private:
Derived decomp_;
// typename Derived::result_type decomp_;
OtherDerived mat_;

public:
Solve(const Derived& decomp, const OtherDerived& mat)
    : decomp_(decomp), mat_(mat) {}
//: decomp_(decomp.get()), mat_(mat) {}

result_type get() const { return decomp_.get().solve(mat_.get()); }
// result_type get() const { return decomp_.solve(mat_.get()); }
};

typedef Holder<Eigen::MatrixXd> MatrixHolder;
typedef Decomp<MatrixHolder> MatrixDecomp;
typedef Solve<MatrixDecomp, MatrixHolder> SimpleSolve;

The following test fails on X.get()

#include "Simple.h"
#include <Eigen/Dense>
#include <iostream>

int main(int, char * []) {
MatrixHolder A(Eigen::MatrixXd::Identity(3, 3));
MatrixHolder B(Eigen::MatrixXd::Random(3, 2));
MatrixDecomp ldlt(A);
SimpleSolve X(ldlt, B);
std::cout << X.get() << std::endl;
return 0;
}

but if you use the commented out lines in the header file, everything works fine. Unfortunately, this moves the evaluation of the decomposition to the construction of the solver, which is not suited to my use. Typically, I want to build a complicated expression expr involving this Solve, and call expr.get() later on.

How can I solve this problem? Is there a general rule to follow, to avoid further related questions?

Answer 1

To avoid useless and expensive copies, the internal solve_retval structure stores the decomposition and right-hand-side by const references. However, the LDLT object created in the Decomp::get function is deleted at the same time this function returns, and thus the solve_retval object references dead objects.

One possible workaround, is to add a Decomp::result_type object in Solve, and initialize it in Solve::get. Moreover, to avoid multiple deep copies, I suggest to use const references for a few attributes as follow:

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

template <class EigenType> class Holder {
public:
  typedef EigenType result_type;

private:
  result_type in_;

public:
  Holder(const EigenType& in) : in_(in) {}
  const result_type& get() const { return in_; }
};

template <class Hold> class Decomp {
public:
  typedef typename Eigen::LDLT
      <typename Eigen::MatrixBase<typename Hold::result_type>::PlainObject>
          result_type;

private:
  const Hold& mat_;
  mutable result_type result_;
  mutable bool init_;

public:
  Decomp(const Hold& mat) : mat_(mat), init_(false) {}

  const result_type& get() const {
    if(!init_) {
      init_ = true;
      result_.compute(mat_.get());
      return result_;
    }
  }
};

template <class Derived, class OtherDerived> class Solve {
public:
  typedef typename Eigen::internal::solve_retval
      <typename Derived::result_type, typename OtherDerived::result_type>
          result_type;

private:
  const Derived& decomp_;
  const OtherDerived& mat_;

public:
  Solve(const Derived& decomp, const OtherDerived& mat)
      : decomp_(decomp), mat_(mat) {}

  result_type get() const {
    return decomp_.get().solve(mat_.get());
  }
};

The general rule is to store heavy objects by const reference (to avoid deep copies) and lightweight expressions by value (to reduce temporary life issues).