Sample code for a variance function

Variance.h

/**
 * @brief Generic class variance which provides a static compute method. This class
 * can work with generic matricies.
 */
template <enum Backend, class Matrix>
struct variance
{
	/** Scalar type */
	typedef typename Matrix::Scalar T;

	/** Calculates unbiased empirical variance estimate given the entries from a matrix. Given a matrix
	 * \f$x\f$ with entries \f$\{x_{11}, ..., x_{mn}\}\f$, this is
	 * \f$\frac{1}{m*n-1}\sum_{i=1}^m\sum_{j=1}^n (x_{ij}-\bar{x})^2\f$ where
	 * \f$\bar x=\frac{1}{mn}\sum_{i=1}^m\sum_{j=1}^n x_{ij}\f$
	 *
	 * @param x matrix of values
	 * @return variance of given values
	 */
	static T compute(Matrix x);
};

/**
 * @brief Specialization of element-wise variance which works with SGMatrix
 * and uses Eigen3 as backend for computing variance.
 */
template <class Matrix>
struct variance<Backend::EIGEN3, Matrix>
{
	/** Scalar type */
	typedef typename Matrix::Scalar T;

	/** Eigen vector type */
	typedef Eigen::Matrix<T,Eigen::Dynamic, Eigen::Dynamic> MatrixXt;

	/** Calculates unbiased empirical variance estimate given the entries from a matrix. Given a matrix
	 * \f$x\f$ with entries \f$\{x_{11}, ..., x_{mn}\}\f$, this is
	 * \f$\frac{1}{m*n-1}\sum_{i=1}^m\sum_{j=1}^n (x_{ij}-\bar{x})^2\f$ where
	 * \f$\bar x=\frac{1}{mn}\sum_{i=1}^m\sum_{j=1}^n x_{ij}\f$
	 *
	 * @param x matrix of values
	 * @return variance of given values
	 */
	static T compute(SGMatrix<T> x)
	{
		REQUIRE(x.num_rows > 0, "Please ensure that m has more than 0 rows.\n")
		REQUIRE(x.num_cols > 0, "Please ensure that m has more than %d columns.\n")

		Eigen::Map<MatrixXt> eigX = x;
		Eigen::Map<MatrixXt> eigSquaredResult = squaredResult;
		SGMatrix<T> squaredResult(x.num_rows, x.num_cols);

		T meanVal = mean<Backend::EIGEN3, SGMatrix<T>>::compute(x);		
		eigSquaredResult.fill(meanVal);
		eigSquaredResult = eigX - eigSquaredResult;
		squaredResult = elementwise_square<Backend::EIGEN3, SGMatrix<T>>::compute(squaredResult);

		return ((T) 1 / (x.num_rows*x.num_cols - 1)) * sum<Backend::EIGEN3, SGMatrix<T>>::compute(squaredResult, false);
	}
};