Persistent expression templates with unique_ptr and matrices

Question

I want to use expression templates to create a tree of objects that persists across statement. Building the tree initially involves some computations with the Eigen linear algebra library. The persistent expression template will have additional methods to compute other quantities by traversing the tree in different ways (but I'm not there yet).

To avoid problems with temporaries going out of scope, subexpression objects are managed through std::unique_ptr. As the expression tree is built, the pointers should be propagated upwards so that holding the pointer for the root object ensures all objects are kept alive. The situation is complicated by the fact that Eigen creates expression templates holding references to temporaries that go out of scope at the end of the statement, so all Eigen expressions must be evaluated while the tree is being constructed.

Below is a scaled-down implementation that seems to work when the val type is an object holding an integer, but with the Matrix type it crashes while constructing the output_xpr object. The reason for the crash seems to be that Eigen's matrix product expression template (Eigen::GeneralProduct) gets corrupted before it is used. However, none of the destructors either of my own expression objects or of GeneralProduct seems to get called before the crash happens, and valgrind doesn't detect any invalid memory accesses.

Any help will be much appreciated! I'd also appreciate comments on my use of move constructors together with static inheritance, maybe the problem is there somewhere.

#include <iostream>
#include <memory>

#include <Eigen/Core>

typedef Eigen::MatrixXi val;

// expression_ptr and derived_ptr: contain unique pointers
// to the actual expression objects

template<class Derived>
struct expression_ptr {
    Derived &&transfer_cast() && {
        return std::move(static_cast<Derived &&>(*this));
    }
};

template<class A>
struct derived_ptr : public expression_ptr<derived_ptr<A>> {
    derived_ptr(std::unique_ptr<A> &&p) : ptr_(std::move(p)) {}
    derived_ptr(derived_ptr<A> &&o) : ptr_(std::move(o.ptr_)) {}

    auto operator()() const {
        return (*ptr_)();
    }

private:
    std::unique_ptr<A> ptr_;
};

// value_xpr, product_xpr and output_xpr: expression templates
// doing the actual work

template<class A>
struct value_xpr {
    value_xpr(const A &v) : value_(v) {}

    const A &operator()() const {
        return value_;
    }

private:
    const A &value_;
};

template<class A,class B>
struct product_xpr {
    product_xpr(expression_ptr<derived_ptr<A>> &&a, expression_ptr<derived_ptr<B>> &&b) :
        a_(std::move(a).transfer_cast()), b_(std::move(b).transfer_cast()) {
    }

    auto operator()() const {
        return a_() * b_();
    }

private:
    derived_ptr<A> a_;
    derived_ptr<B> b_;
};

// Top-level expression with a matrix to hold the completely
// evaluated output of the Eigen calculations
template<class A>
struct output_xpr {
    output_xpr(expression_ptr<derived_ptr<A>> &&a) :
        a_(std::move(a).transfer_cast()), result_(a_()) {}

    const val &operator()() const {
        return result_;
    }

private:
    derived_ptr<A> a_;
    val result_;
};

// helper functions to create the expressions

template<class A>
derived_ptr<value_xpr<A>> input(const A &a) {
    return derived_ptr<value_xpr<A>>(std::make_unique<value_xpr<A>>(a));
}

template<class A,class B>
derived_ptr<product_xpr<A,B>> operator*(expression_ptr<derived_ptr<A>> &&a, expression_ptr<derived_ptr<B>> &&b) {
    return derived_ptr<product_xpr<A,B>>(std::make_unique<product_xpr<A,B>>(std::move(a).transfer_cast(), std::move(b).transfer_cast()));
}

template<class A>
derived_ptr<output_xpr<A>> eval(expression_ptr<derived_ptr<A>> &&a) {
    return derived_ptr<output_xpr<A>>(std::make_unique<output_xpr<A>>(std::move(a).transfer_cast()));
}

int main() {
    Eigen::MatrixXi mat(2, 2);
    mat << 1, 1, 0, 1;
    val one(mat), two(mat);
    auto xpr = eval(input(one) * input(two));
    std::cout << xpr() << std::endl;
    return 0;
}

Answer 1

Your problem appears to be that you are using someone else's expression templates, and storing the result in an auto.

(This happens in product_xpr<A>::operator(), where you call *, which if I read it right, is an Eigen multiplication that uses expression templates).

Expression templates are often designed to presume the entire expression will occur on a single line, and it will end with a sink type (like a matrix) that causes the expression template to be evaluated.

In your case, you have a*b expression template, which is then used to construct an expression template return value, which you later evaluate. The lifetime of temporaries passed to * in a*b are going to be over by the time you reach the sink type (matrix), which violates what the expression templates expect.

I am struggling to come up with a solution to ensure that all temporary objects have their lifetime extended. One thought I had was some kind of continuation passing style, where instead of calling:

Matrix m = (a*b);

you do

auto x = { do (a*b) pass that to (cast to matrix) }

replace

auto operator()() const {
    return a_() * b_();
}

with

template<class F>
auto operator()(F&& f) const {
    return std::forward<F>(f)(a_() * b_());
}

where the "next step' is passed to each sub-expression. This gets trickier with binary expressions, in that you have to ensure that the evaluation of the first expression calls code that causes the second sub expression to be evaluated, and then the two expressions are combined, all in the same long recursive call stack.

I am not proficient enough in continuation passing style to untangle this knot completely, but it is somewhat popular in the functional programming world.

Another approach would be to flatten your tree into a tuple of optionals, then construct each optional in the tree using a fancy operator(), and manually hook up the arguments that way. Basically do manual memory management of the intermediate values. This will work if the Eigen expression templates are either move-aware or do not have any self-pointers, so that moving at the point of construction doesn't break things. Writing that would be challenging.

Answer 2

Continuation passing style, suggested by Yakk, solves the problem and isn't too insane (not more insane than template metaprogramming in general anyhow). The double lambda evaluation for the arguments of binary expressions can be tucked away in a helper function, see binary_cont in the code below. For reference, and since it's not entirely trivial, I'm posting the fixed code here.

If somebody understands why I had to put a const qualifier on the F type in binary_cont, please let me know.

#include <iostream>
#include <memory>

#include <Eigen/Core>

typedef Eigen::MatrixXi val;

// expression_ptr and derived_ptr: contain unique pointers
// to the actual expression objects

template<class Derived>
struct expression_ptr {
    Derived &&transfer_cast() && {
        return std::move(static_cast<Derived &&>(*this));
    }
};

template<class A>
struct derived_ptr : public expression_ptr<derived_ptr<A>> {
    derived_ptr(std::unique_ptr<A> &&p) : ptr_(std::move(p)) {}
    derived_ptr(derived_ptr<A> &&o) = default;

    auto operator()() const {
        return (*ptr_)();
    }

    template<class F>
    auto operator()(F &&f) const {
        return (*ptr_)(std::forward<F>(f));
    }

private:
    std::unique_ptr<A> ptr_;
};

template<class A,class B,class F>
auto binary_cont(const derived_ptr<A> &a_, const derived_ptr<B> &b_, const F &&f) {
    return a_([&b_, f = std::forward<const F>(f)] (auto &&a) {
        return b_([a = std::forward<decltype(a)>(a), f = std::forward<const F>(f)] (auto &&b) {
            return std::forward<const F>(f)(std::forward<decltype(a)>(a), std::forward<decltype(b)>(b));
        });
    });
}

// value_xpr, product_xpr and output_xpr: expression templates
// doing the actual work

template<class A>
struct value_xpr {
    value_xpr(const A &v) : value_(v) {}

    template<class F>
    auto operator()(F &&f) const {
        return std::forward<F>(f)(value_);
    }

private:
    const A &value_;
};

template<class A,class B>
struct product_xpr {
    product_xpr(expression_ptr<derived_ptr<A>> &&a, expression_ptr<derived_ptr<B>> &&b) :
        a_(std::move(a).transfer_cast()), b_(std::move(b).transfer_cast()) {
    }

    template<class F>
    auto operator()(F &&f) const {
        return binary_cont(a_, b_,
            [f = std::forward<F>(f)] (auto &&a, auto &&b) {
                return f(std::forward<decltype(a)>(a) * std::forward<decltype(b)>(b));
            });
    }

private:
    derived_ptr<A> a_;
    derived_ptr<B> b_;
};

template<class A>
struct output_xpr {
    output_xpr(expression_ptr<derived_ptr<A>> &&a) :
            a_(std::move(a).transfer_cast()) {
        a_([this] (auto &&x) { this->result_ = x; });
    }

    const val &operator()() const {
        return result_;
    }

private:
    derived_ptr<A> a_;
    val result_;
};

// helper functions to create the expressions

template<class A>
derived_ptr<value_xpr<A>> input(const A &a) {
    return derived_ptr<value_xpr<A>>(std::make_unique<value_xpr<A>>(a));
}

template<class A,class B>
derived_ptr<product_xpr<A,B>> operator*(expression_ptr<derived_ptr<A>> &&a, expression_ptr<derived_ptr<B>> &&b) {
    return derived_ptr<product_xpr<A,B>>(std::make_unique<product_xpr<A,B>>(std::move(a).transfer_cast(), std::move(b).transfer_cast()));
}

template<class A>
derived_ptr<output_xpr<A>> eval(expression_ptr<derived_ptr<A>> &&a) {
    return derived_ptr<output_xpr<A>>(std::make_unique<output_xpr<A>>(std::move(a).transfer_cast()));
}

int main() {
    Eigen::MatrixXi mat(2, 2);
    mat << 1, 1, 0, 1;
    val one(mat), two(mat), three(mat);
    auto xpr = eval(input(one) * input(two) * input(one) * input(two));
    std::cout << xpr() << std::endl;
    return 0;
}