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 L(3,3), P(3,3);
Eigen::Matrix L(3,K), P(3,K);
Eigen::Matrix L(N,K), P(N,K);