关于数学库Eigen两个常用函数实现（伪逆矩阵和协方差矩阵）

//求解协方差矩阵
Eigen::MatrixXd& cov(Eigen::MatrixXd& outMatrix, const Eigen::MatrixXd& inMatrix)
{
	Eigen::MatrixXd meanVec = inMatrix.colwise().mean();
	Eigen::RowVectorXd meanVecRow(Eigen::RowVectorXd::Map(meanVec.data(), inMatrix.cols()));
	Eigen::MatrixXd zeroMeanMat = inMatrix;
	zeroMeanMat.rowwise() -= meanVecRow;
	if (inMatrix.rows() == 1)
	{
		outMatrix = (zeroMeanMat.adjoint() * zeroMeanMat) / double(inMatrix.rows());
	}
	else
	{
		outMatrix = (zeroMeanMat.adjoint() * zeroMeanMat) / double(inMatrix.rows() - 1);
	}
	return outMatrix;
}
////////////////////////////////////////////////////////////////////////////////////////////////////
//伪逆矩阵(Moore-Penrose pseudoinverse)A定义：
//A+=VD+UT,其中，U，D和V是矩阵A奇异值分解后得到的矩阵。对角矩阵D的伪逆D+是非零元素取倒数之后再转置得到的。
//
Eigen::MatrixXd& pinv(Eigen::MatrixXd& outMatrix, const Eigen::MatrixXd& inMatrix)
{
	double pinvtoler = 1.e-6; // choose your tolerance wisely!
	Eigen::JacobiSVD<Eigen::MatrixXd> svd(inMatrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
	Eigen::VectorXd singularValues_inv = svd.singularValues();
	Eigen::VectorXd sv = svd.singularValues();
	for (Eigen::Index i = 0; i < svd.cols(); ++i)
	{
		if (sv(i) > pinvtoler)
		{
			singularValues_inv(i) = 1.0 / sv(i);
		}
		else
		{
			singularValues_inv(i) = 0;
		}
	}
	outMatrix = (svd.matrixV() * singularValues_inv.asDiagonal() * svd.matrixU().transpose());
	return outMatrix;
}