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

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);