协方差矩阵介绍及C++/OpenCV/Eigen的三种实现

函数f(x)关于某分布P(x)的期望(expectation)或者期望值(expected value)是指，当x由P产生，f作用于x时，f(x)的平均值。对于离散型随机变量，这可以通过求和得到：


对于连续型随机变量可以通过求积分得到：

当概率分布在上下文中指明时，我们可以只写出期望作用的随机变量的名称来进行简化，例如Ex[f(x)]。如果期望作用的随机变量也很明确，我们可以完全不写脚标，就像E[f(x)]。默认地，我们假设E[·]表示对方括号内的所有随机变量的值求平均。期望是线性的。 方差(variance)衡量的是当我们对x依据它的概率分布进行采样时，随机变量x的函数值会呈现多大的差异：Var(f(x))=E[(f(x)-E[f(x)])2] 当方差很小时，f(x)的值形成的簇比较接近它们的期望值。方差的平方根被称为标准差(standard deviation)。 协方差(covariance)在某种意义上给出了两个变量线性相关性的强度以及这些变量的尺度：Cov(f(x),g(y))=E[(f(x)-E[f(x)])(g(y)-E[g(y)])] 协方差的绝对值如果很大则意味着变量值变化很大并且它们同时距离各自的均值很远。如果协方差是正的，那么两个变量都倾向于同时取得相对较大的值。如果协方差是负的，那么其中一个变量倾向于取得相对较大的值的同时，另一个变量倾向于取得相对较小的值，反之亦然。其它的衡量指标如相关系数(correlation)将每个变量的贡献归一化，为了只衡量变量的相关性而不受各个变量尺度大小的影响。 协方差和相关性是有联系的，但实际上不同的概念。它们是有联系的，因为两个变量如果相互独立那么它们的协方差为零，如果两个变量的协方差不为零那么它们一定是相关的。然而，独立性又和协方差完全不同的性质。两个变量如果协方差为零，它们之间一定没有线性关系。独立性是比零协方差的要求更强，因为独立性还排除了非线性的关系。两个变量相互依赖但是具有零协方差是可能的。
随机向量x∈Rn的协方差矩阵(covariance matrix)是一个n*n的矩阵，并且满足：Cov(x)i,j=Cov(xi,xj) 协方差矩阵的对角元是方差：Cov(xi,xi)=Var(xi). 期望值：在概率论和统计学中，一个离散性随机变量的期望值(或数学期望、或均值，亦简称期望)是试验中每次可能的结果乘以其结果概率的总和。换句话说，期望值像是随机试验在同样的机会下重复多次，所有那些可能状态平均的结果，便基本上等同”期望值”所期望的数。需要注意的是，期望值并不一定等同于常识中的”期望”------“期望值”也许与每一个结果都不相等。(换句话说，期望值是该变量输出值的平均数。期望值并不一定包含于变量的输出值集合里。) 期望值数学定义：如果X是在概率区间(Ω,P)中的随机变量，那么它的期望值E[X]的定义是：E[X]=∫ΩXdP，并不是每一个随机变量都有期望值的，因为有的时候这个积分不存在。如果两个随机变量的分布相同，那它们的期望值也相同。 如果X是离散的随机变量，输出值为x1,x2,…,和输出值相应的概率为p1,p2,…(概率和为1)。若级数∑ipixi绝对收敛，那么期望值E[X]是一个无限数列的和：E[X]=∑pixi 期望值性质： (1)、期望值E是线性函数：E[aX+bY]=aE[X]+bE[Y]，X和Y为在同一概率空间的两个随机变量(可以独立或者非独立),a和b为任意实数。 (2)、一般的说，一个随机变量的函数的期望值并不等于这个随机变量的期望值的函数。 (3)、在一般情况下，两个随机变量的积的期望值不等于两个随机变量的期望值的积。特殊情况是当这两个随机变量是相互独立的时候E[XY]=E[X]E[Y](也就是说一个随机变量的输出不会影响另一个随机变量的输出)。 在统计学中，当估算变量的期望值时，经常用到的方法是重复测量此变量的值，然后用所得数据的平均值来作为此变量的期望值的估计。在概率分布中，期望值和方差或标准差是一种分布的重要特征。 期望值也可以通过方差计算公式来计算方差：Var(X)=E[X2]-(E[X])2。 方差(Variance)：在概率论和统计学中，一个随机变量的方差描述的是它的离散程度，也就是该变量离其期望值的距离。一个实随机变量的方差也称为它的二阶炬或二阶中心动差，恰巧也是它的二阶累积量。这里把复杂说白了，就是将各个误差将之平方(而非取绝对值)，使之肯定为正数，相加之后再除以总数，透过这样的方式来算出各个数据分布、零散(相对中心点)的程度。继续延伸的话，方差的算术平方根称为该随机变量的标准差(此为相对各个数据点间)。 方差定义：设X为服从分布F的随机变量，如果E[X]是随机变数X的期望值(平均数μ=E[X])，随机变量X或者分布F的方差为: Var(X)=E[(X-μ)2] ,这个定义涵盖了连续、离散、或两者都有的随机变数。方差亦可当作是随机变数与自己本身的协方差：Var(X)=Cov(X,X)。 方差表示式展开成为：Var(X)=E[X2-2XE[X]+(E[X])2]=E[X2]-2E[X]E[X]+(E[X])2=E[X2]-(E[X])2,上述的表示式可记为”平方的平均减掉平均的平方”。 方差不会是负的。 标准差(Standard Deviation, SD)，数学符号σ(sigma)，在概率统计中最常使用作为测量一组数值的离散程度之用。标准差定义：为方差开算术平方根，反映组内个体间的离散程度。标准差与期望值之比为标准离差率。 简单来说，标准差是一组数值自平均值分散开来的程度的一种测量观念。一个较大的标准差，代表大部分的数值和其平均值之间差异较大；一个较小的标准差，代表这些数值较接近平均值。 总体的标准差：基本定义：

协方差(Covariance)：在概率论和统计学中用于衡量两个变量的总体误差。而方差是协方差的一种特殊情况，即当两个变量是相同的情况。期望值分别为E(X)=μ与E(Y)=ν的两个实数随机变量X与Y之间的协方差定义为：cov(X,Y)=E((X-μ)(Y-ν)),cov(X,Y)=E(X·Y)- μν. 协方差表示的是两个变量的总体的误差，这与只表示一个变量误差的方差不同。如果两个变量的变化趋势一致，也就是说如果其中一个大于自身的期望值，另外一个也大于自身的期望值，那么两个变量之间的协方差就是正值。如果两个变量的变化趋势相反，即其中一个大于自身的期望值，另外一个却小于自身的期望值，那么两个变量之间的协方差就是负值。 如果X与Y是统计独立的，那么二者之间的协方差就是0，这是因为E(X·Y)=E(X)·E(Y)= μν. 但是反过来并不成立，即如果X与Y的协方差为0，二者并不一定是统计独立的。取决于协方差的相关性η：

更准确地说是线性相关性，是一个衡量线性独立的无量刚数，其取值在[-1,+1]之间。相关性η=1时称为”完全线性相关”(相关性η=-1时称为”完全线性负相关”)，此时将Yi对Xi作Y-X散点图，将得到一组精确排列在直线上的点；相关性数值介于-1到1之间时，其绝对值越接近1表明线性相关性越好，作散点图得到的点的排布越接近一条直线。 相关性为0(因而协方差也为0)的两个随机变量又被称为是不相关的，或者更准确地说叫做”线性无关”、”线性不相关”，这仅仅表明X与Y两随机变量之间没有线性相关性，并非表示它们之间一定没有任何内在的(非线性)函数关系，和前面所说的”X、Y二者并不一定是统计独立的”说法一致。 如果X与Y是实数随机变量，a与b不是随机变量，那么根据协方差的定义可以得到：cov(X,X)=var(X), cov(X,Y)=cov(Y,X),cov(aX,bY)=ab cov(X,Y), 对于随机变量序列X1,…,Xn与Y1,…,Ym，有：

对于随机变量序列X1,…,Xn，有：

两个向量变量的协方差cov(X,Y)与cov(Y,X)互为转置矩阵。 协方差矩阵：在统计学与概率论中，协方差矩阵(也称离差矩阵、方差-协方差矩阵)是一个矩阵，其i,j位置的元素是第i个与第j个随机向量(即随机变量构成的向量)之间的协方差。这是从标量随机变量到高维度随机向量的自然推广。 假设X是以n个随机变数(其中的每个随机变数也是一个向量，是一个列向量)组成的行向量，X=[X1,X2,…,Xn].并且μi是其第i个列向量Xi中所有元素的平均值，即μi=E(Xi)。协方差矩阵的第i,j项(第i,j项是一个协方差)被定义为如下形式：∑ij=cov(Xi,Yj)=E[(Xi-μi)(Xj-μj)T] 而协方差矩阵为：

矩阵中的第(i,j)个元素是Xi与Xj的协方差。这个概念是对于标量随机变数方差的一般化推广。 协方差矩阵性质：∑=E[(X-E[X])(X-E[X])T]与μ=E(X)满足下边的基本性质： (1)、∑=E(XXT)- μμT (2)、∑是半正定的和对称的矩阵 (3)、var(aTX)=aTvar(X)a (4)、∑≥0 (5)、var(AX+a)=Avar(X)AT (6)、cov(X,Y)=cov(Y,X)T (7)、cov(X1+X2,Y)=cov(X1,Y)+cov(X2,Y) (8)、若p=q，则有cov(X+Y)=var(X)+cov(X,Y)+cov(Y,X)+var(Y) (9)、cov(AX,BX)=Acov(X,X)BT (10)、若X与Y是独立的，则有cov(X,Y)=0 (11)、∑=∑T
其中X,X1与X2是随机(p*1)向量，Y是随机(q*1)向量，a是(p*1)向量，A与B是(q*p)矩阵。
以上内容摘自：《深度学习中文版》和 维基百科
以下是分别采用C++和OpenCV实现的计算协方差矩阵code：
#include "funset.hpp"
#include <math.h>
#include <iostream>
#include <string>
#include <vector>
#include <opencv2/opencv.hpp>
#include "common.hpp"

// =============================== 计算协方差矩阵 ============================
// 按行存储，以行为向量,输入矩阵mat为m行n列，则协方差矩阵covar为n行n列对称矩阵，均值mean为1行n列
template<typename _Tp>
int calcCovarMatrix(const std::vector<std::vector<_Tp>>& mat, std::vector<std::vector<_Tp>>& covar, std::vector<_Tp>& mean, bool scale = false)
{
const int rows = mat.size();
const int cols = mat[0].size();
const int nsamples = rows;
double scale_ = 1.;
if (scale) scale_ = 1. / (nsamples /*- 1*/);

covar.resize(cols);
for (int i = 0; i < cols; ++i)
covar[i].resize(cols, (_Tp)0);
mean.resize(cols, (_Tp)0);

for (int w = 0; w < cols; ++w) {
for (int h = 0; h < rows; ++h) {
mean[w] += mat[h][w];
}
}

for (auto& value : mean) {
value = 1. / rows * value;
}

for (int i = 0; i < cols; ++i) {
std::vector<_Tp> col_buf(rows, (_Tp)0);
for (int k = 0; k < rows; ++k)
col_buf[k] = mat[k][i] - mean[i];

for (int j = 0; j < cols; ++j) {
double s0 = 0;
for (int k = 0; k < rows; ++k) {
s0 += col_buf[k] * (mat[k][j] - mean[j]);
}
covar[i][j] = (_Tp)(s0 * scale_);
}
}

return 0;
}

int test_calcCovarMatrix()
{
std::vector<std::vector<float>> vec{ { 1.2f, 2.5f, 5.6f, -2.5f },
{ -3.6f, 9.2f, 0.5f, 7.2f },
{ 4.3f, 1.3f, 9.4f, -3.4f } };
const int rows{ 3 }, cols{ 4 };

fprintf(stderr, "source matrix:\n");
fbc::print_matrix(vec);

fprintf(stderr, "\nc++ implement calculate covariance matrix:\n");
std::vector<std::vector<float>> covar1;
std::vector<float> mean1;
if (fbc::calcCovarMatrix(vec, covar1, mean1, false/*true*/) != 0) {
fprintf(stderr, "C++ implement calcCovarMatrix fail\n");
return -1;
}

fprintf(stderr, "print covariance matrix: \n");
fbc::print_matrix(covar1);
fprintf(stderr, "print mean: \n");
fbc::print_matrix(mean1);

fprintf(stderr, "\nc++ opencv calculate covariance matrix:\n");
cv::Mat mat(rows, cols, CV_32FC1);
for (int y = 0; y < rows; ++y) {
for (int x = 0; x < cols; ++x) {
mat.at<float>(y, x) = vec.at(y).at(x);
}
}
//std::cout << mat << std::endl;

//std::cout << "mat:" << std::endl << mat << std::endl;
//cv::Mat means(1, mat.cols, mat.type(), cv::Scalar::all(0));
//for (int i = 0; i < mat.cols; i++)
// means.col(i) = (cv::sum(mat.col(i)) / mat.rows).val[0];
//std::cout << "means:" << std::endl << means << std::endl;
//cv::Mat tmp = cv::repeat(means, mat.rows, 1);
//cv::Mat mat2 = mat - tmp;
//cv::Mat covar = (mat2.t()*mat2) / (mat2.rows /*- 1*/); // （X'*X)/n-1
//std::cout << "covar:" << std::endl << covar << std::endl;

cv::Mat covar2, mean2;
cv::calcCovarMatrix(mat, covar2, mean2, CV_COVAR_NORMAL | CV_COVAR_ROWS/* | CV_COVAR_SCALE*/, CV_32FC1);
fprintf(stderr, "print covariance matrix: \n");
fbc::print_matrix(covar2);
fprintf(stderr, "print mean: \n");
fbc::print_matrix(mean2);

return 0;
} 执行结果如下：


以下是采用Eigen实现的计算协方差矩阵code：

#include "funset.hpp"
#include <math.h>
#include <iostream>
#include <vector>
#include <string>
#include <opencv2/opencv.hpp>
#include <Eigen/Dense>
#include "common.hpp"

int test_calcCovarMatrix()
{
// reference: https://stackoverflow.com/questions/15138634/eigen-is-there-an-inbuilt-way-to-calculate-sample-covariance std::vector<std::vector<float>> vec{ { 1.2f, 2.5f, 5.6f, -2.5f },
{ -3.6f, 9.2f, 0.5f, 7.2f },
{ 4.3f, 1.3f, 9.4f, -3.4f } };
const int rows{ 3 }, cols{ 4 };

std::vector<float> vec_;
for (int i = 0; i < rows; ++i) {
vec_.insert(vec_.begin() + i * cols, vec[i].begin(), vec[i].end());
}
Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> m(vec_.data(), rows, cols);

fprintf(stderr, "source matrix:\n");
std::cout << m << std::endl;

fprintf(stdout, "\nEigen implement:\n");
const int nsamples = rows;
float scale = 1. / (nsamples /*- 1*/);

Eigen::MatrixXf mean = m.colwise().mean();
std::cout << "print mean: " << std::endl << mean << std::endl;

Eigen::MatrixXf tmp(rows, cols);
for (int y = 0; y < rows; ++y) {
for (int x = 0; x < cols; ++x) {
tmp(y, x) = m(y, x) - mean(0, x);
}
}
//std::cout << "tmp: " << std::endl << tmp << std::endl;

Eigen::MatrixXf covar = (tmp.adjoint() * tmp) /*/ float(nsamples - 1)*/;
std::cout << "print covariance matrix: " << std::endl << covar << std::endl;

return 0;
}

执行结果如下：


由以上结果可见：C++、OpenCV、Eigen实现结果是一致的。 注：以上是按照行存储，且无scale时的结果，如果要进行scale，在OpenCV中实现时是除以行数(nsamples)，而在其它公式中介绍是除以行数-1。
GitHub： https://github.com/fengbingchun/NN_Test https://github.com/fengbingchun/Eigen_Test