Assuming your samples are stored row-wise in a matrix D, then you can do:
VectorXd N = D.rowwise().squaredNorm();
MatrixXd S = N.replicate(1,n) + N.transpose().replicate(n,1);
S.noalias() -= 2. * D * D.transpose();
S = S.array().sqrt();
This exploits the fact that |x-y|²=x²+y²-2x'y
. The noalias() statement is just an optimization to Eigen there is no risk of aliasing in this product, thus no temporary is needed. The .array() statement switches to the array world where all functions are applied coefficient-wise.