Probably won't ever catch Fortran, but using Eigen templates for reductions really opens the door for compile time optimizations; e.g. these are all reductions that do the same thing

    const float residual = (L.array() * P.array()).colwise().sum().square().mean();
    const float residual = L.cwiseProduct(P).array().colwise().sum().array().square().mean();
    const float residual = (L.transpose() * P).diagonal().array().square().mean();

The compiler can optimize using static information, e.g. these would all be handled differently for the following types

    Eigen::Matrix<float,3,3> L(3,3), P(3,3);
    Eigen::Matrix<float,3,Eigen::Dynamic> L(3,K), P(3,K);
    Eigen::Matrix<float,Eigen::Dynamic,Eigen::Dynamic> L(N,K), P(N,K);