Eigen deduce Template Expression type

Question

In short: I am trying to deduce the return value (which is an expression template) of a function that performs a unary operation on a matrix.

In this case the operation is to compute a covariance matrix.

I followed the Eigen documentation here: https://eigen.tuxfamily.org/dox/TopicCustomizing_NullaryExpr.html and created a minimal example of multiplying a matrix by two and returning its result. Code snippet below shows a working example. The key point for me is that the expression is not evaluated to an intermediate result so I do not want to return something like Eigen::MatrixBase<Derived>.

template<typename ArgType>
struct times_two_helper {
    template<typename Derived>
    static auto TimesTwo(Eigen::MatrixBase<Derived> const& mat) {
        return mat + mat;
    }
    using ResultType = Eigen::Matrix<typename ArgType::Scalar,
            ArgType::ColsAtCompileTime,
            ArgType::RowsAtCompileTime,
            ArgType::Options,
            ArgType::MaxColsAtCompileTime,
            ArgType::MaxRowsAtCompileTime>;
    using ExpressionType = decltype(TimesTwo(std::declval<ArgType>()));
};
template<typename ArgType>
struct times_two_functor {
    using ResultType = typename times_two_helper<ArgType>::ResultType;
    using ExpressionType = typename times_two_helper<ArgType>::ExpressionType;
    ArgType const& arg_;

    // Here I want to store the expression without evaluating it!
    ExpressionType expression_;

public:
    times_two_functor(ArgType const& arg)
        : arg_{arg}
        , expression_{times_two_helper<ArgType>::TimesTwo(arg)}
    {}
    typename ArgType::Scalar operator() (Eigen::Index row, Eigen::Index col) const {
        return expression_(row, col);
    }
};
template <class ArgType>
Eigen::CwiseNullaryOp<times_two_functor<ArgType>, typename times_two_helper<ArgType>::ResultType>
TimesTwo(Eigen::MatrixBase<ArgType> const& arg)         {
    using ResultType = typename times_two_helper<ArgType>::ResultType;
    return ResultType::NullaryExpr(arg.rows(), arg.cols(), times_two_functor<ArgType>(arg.derived()));
}

Used like this:

TEST(Stat, TimesTwo) {
    Eigen::Matrix<double, 3, 3> input;
    input << 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0;

    Eigen::Matrix<double, 3, 3> result = TimesTwo(input);
    std::cout << result << "\n";
}

But when I try to do the same thing for a covariance matrix it deduces the wrong type for ExpressionType (error message is given at the bottom).

template<class ArgType>
struct covariance_helper {
    template<typename Derived>
    static auto Covariance(Eigen::MatrixBase<Derived> const& mat) {
        auto centered = mat.rowwise() - mat.colwise().mean();
        return (centered.adjoint() * centered) / double(mat.rows() - 1);
    }

    using ResultType = Eigen::Matrix<typename ArgType::Scalar,
            ArgType::ColsAtCompileTime,
            ArgType::ColsAtCompileTime,
            ArgType::Options,
            ArgType::MaxColsAtCompileTime,
            ArgType::MaxColsAtCompileTime>;
    using ExpressionType = decltype(covariance_helper::Covariance(std::declval<ArgType>()));
};
template<class ArgType>
class covariance_functor {
    using ResultType = typename times_two_helper<ArgType>::ResultType;
    using ExpressionType = typename times_two_helper<ArgType>::ExpressionType;
    const ArgType &mat_;
    ExpressionType expression_;
public:
    covariance_functor(const ArgType& arg)
        : mat_{arg}
        , expression_(covariance_helper<ArgType>::Covariance(arg))
    {}
    typename ArgType::Scalar operator() (Eigen::Index row, Eigen::Index col) const {
        return expression_(row, col);
    }
};
template <class ArgType>
Eigen::CwiseNullaryOp<covariance_functor<ArgType>, typename covariance_helper<ArgType>::ResultType>
Covariance(Eigen::MatrixBase<ArgType> const& arg) {
    using ResultType = typename covariance_helper<ArgType>::ResultType;
    return ResultType::NullaryExpr(arg.cols(), arg.cols(), covariance_functor<ArgType>(arg.derived()));
}

Calling it like this:

TEST(Stat, Covariance) {
    constexpr double kPrecision = 1e-12;
    // source: https://www.itl.nist.gov/div898/handbook/pmc/section5/pmc541.htm
    Eigen::Matrix<double, 5, 3> mat;
    mat << 4.0, 2.0, 0.6, 4.2, 2.1, 0.59, 3.9, 2.0, 0.58, 4.3, 2.1, 0.62, 4.1, 2.2, 0.63;

    Eigen::Matrix<double, 3, 3> cov_expected;
    cov_expected << 0.025, 0.0075, 0.00175, 0.0075, 0.0070, 0.00135, 0.00175, 0.00135, 0.00043;

    Covariance(mat);
}

gives me the following error message:

In file included from /test/test_stat.cpp:8:0:
/lib/core/stat.hpp: In instantiation of ‘covariance_functor<ArgType>::covariance_functor(const ArgType&) [with ArgType = Eigen::Matrix<double, 5, 3>]’:
/lib/core/stat.hpp:138:60:   required from ‘Eigen::CwiseNullaryOp<covariance_functor<ArgType>, typename covariance_helper<ArgType>::ResultType> Covariance(const Eigen::MatrixBase<Derived>&) [with ArgType = Eigen::Matrix<double, 5, 3>; typename covariance_helper<ArgType>::ResultType = Eigen::Matrix<double, 3, 3>]’
/test/test_stat.cpp:19:29:   required from here
/lib/core/stat.hpp:127:66: error: no matching function for call to ‘Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<double, double>, const Eigen::Matrix<double, 5, 3>, const Eigen::Matrix<double, 5, 3> >::CwiseBinaryOp(Eigen::CwiseBinaryOp<Eigen::internal::scalar_quotient_op<double, double>, const Eigen::Product<Eigen::Transpose<const Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, const Eigen::Matrix<double, 5, 3>, const Eigen::Replicate<Eigen::PartialReduxExpr<const Eigen::Matrix<double, 5, 3>, Eigen::internal::member_mean<double>, 0>, 5, 1> > >, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, const Eigen::Matrix<double, 5, 3>, const Eigen::Replicate<Eigen::PartialReduxExpr<const Eigen::Matrix<double, 5, 3>, Eigen::internal::member_mean<double>, 0>, 5, 1> >, 0>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, const Eigen::Matrix<double, 3, 3> > >)’
         , expression_(covariance_helper<ArgType>::Covariance(arg))

What am I missing here? Is the thing I am trying to do possible? Is it reasonable to want this in the context of delayed evaluation (specific to covariance matrices)?

This is an exercise for me (not homework) so I really want to know if it is possible to deduce the full type of an expression template in this way.

I am of course also interested in the practicality of this approach but to a lesser extend.

EDIT:

Since this might help someone, here is a working version.

template<class ArgType>
struct covariance_helper {
    template<typename Derived>
    static auto Covariance(Eigen::MatrixBase<Derived> const& mat) {
        auto centered = mat.rowwise() - mat.colwise().mean();
        return (centered.adjoint() * centered) / double(mat.rows() - 1);
    }

    using ResultType = Eigen::Matrix<typename ArgType::Scalar,
            ArgType::ColsAtCompileTime,
            ArgType::ColsAtCompileTime,
            ArgType::Options,
            ArgType::MaxColsAtCompileTime,
            ArgType::MaxColsAtCompileTime>;
    using ExpressionType = decltype(covariance_helper::Covariance(std::declval<Eigen::MatrixBase<ArgType> const&>()));
};
template<class ArgType>
class covariance_functor {
    using ResultType = typename covariance_helper<ArgType>::ResultType;
    using ExpressionType = typename covariance_helper<ArgType>::ExpressionType;
    const ArgType &mat_;
    ExpressionType expression_;
public:
    covariance_functor(const ArgType& arg)
        : mat_{arg}
        , expression_(covariance_helper<ArgType>::Covariance(arg))
    {}
    typename ArgType::Scalar operator() (Eigen::Index row, Eigen::Index col) const {
        return expression_(row, col);
    }
};
template <class ArgType>
Eigen::CwiseNullaryOp<covariance_functor<ArgType>, typename covariance_helper<ArgType>::ResultType>
Covariance(Eigen::MatrixBase<ArgType> const& arg) {
    using ResultType = typename covariance_helper<ArgType>::ResultType;
    return ResultType::NullaryExpr(arg.cols(), arg.cols(), covariance_functor<ArgType>(arg.derived()));
}

Call it like this:

TEST(Stat, Covariance) {
    constexpr double kPrecision = 1e-12;
    // source: https://www.itl.nist.gov/div898/handbook/pmc/section5/pmc541.htm
    Eigen::Matrix<double, 5, 3> mat;
    mat << 4.0, 2.0, 0.6, 4.2, 2.1, 0.59, 3.9, 2.0, 0.58, 4.3, 2.1, 0.62, 4.1, 2.2, 0.63;

    Eigen::Matrix<double, 3, 3> cov_expected;
    cov_expected << 0.025, 0.0075, 0.00175, 0.0075, 0.0070, 0.00135, 0.00175, 0.00135, 0.00043;

    Eigen::MatrixXd cov = Covariance(mat);
    ASSERT_NEAR((cov - cov_expected).norm(), 0.0, kPrecision);
}

Or a simpler one:

template<typename Derived>
struct Covariance {

    using MBase = Eigen::MatrixBase<Derived> const&;

    constexpr static auto compute(MBase mat) {
        auto centered = mat.rowwise() - mat.colwise().mean();
        return (centered.adjoint() * centered) / double(mat.rows() - 1);
    }

    using ResultExpr = decltype(compute(std::declval<MBase>()));

    Covariance(MBase mat)
            : result_expr{compute(mat)}
    {}

    typename Derived::Scalar operator() (Eigen::Index row, Eigen::Index col) const {
        return result_expr(row, col);
    }

    ResultExpr result_expr;
};
TEST(Stat, EigenUnaryExpr) {
    constexpr double kPrecision = 1e-12;
    // source: https://www.itl.nist.gov/div898/handbook/pmc/section5/pmc541.htm
    Eigen::Matrix<double, 5, 3> mat;
    mat << 4.0, 2.0, 0.6, 4.2, 2.1, 0.59, 3.9, 2.0, 0.58, 4.3, 2.1, 0.62, 4.1, 2.2, 0.63;

    Eigen::Matrix<double, 3, 3> cov_expected;
    cov_expected << 0.025, 0.0075, 0.00175, 0.0075, 0.0070, 0.00135, 0.00175, 0.00135, 0.00043;

    Covariance cov{mat};
    using ResultExpr = decltype(cov)::ResultExpr;
    auto cov_expr = ResultExpr::NullaryExpr(cov.result_expr.rows(), cov.result_expr.cols(), cov);

    std::cerr << cov_expr << "\n";
    ASSERT_NEAR((cov_expr - cov_expected).norm(), 0.0, kPrecision);
}