协方差的意义和计算公式整理

学过概率统计的孩子都知道，统计里最基本的概念就是样本的均值，方差，或者再加个标准差。首先我们给你一个含有n个样本的集合，依次给出这些概念的公式描述，这些高中学过数学的孩子都应该知道吧，一带而过。

很显然，均值描述的是样本集合的中间点，它告诉我们的信息是很有限的，而标准差给我们描述的则是样本集合的各个样本点到均值的距离之平均。以这两个集合为例，[0，8，12，20]和[8，9，11，12]，两个集合的均值都是10，但显然两个集合差别是很大的，计算两者的标准差，前者是8.3，后者是1.8，显然后者较为集中，故其标准差小一些，标准差描述的就是这种“散布度”。之所以除以n-1而不是除以n，是因为这样能使我们以较小的样本集更好的逼近总体的标准差，即统计上所谓的“无偏估计”。而方差则仅仅是标准差的平方。

为什么需要协方差？

上面几个统计量看似已经描述的差不多了，但我们应该注意到，标准差和方差一般是用来描述一维数据的，但现实生活我们常常遇到含有多维数据的数据集，最简单的大家上学时免不了要统计多个学科的考试成绩。面对这样的数据集，我们当然可以按照每一维独立的计算其方差，但是通常我们还想了解更多，比如，一个男孩子的猥琐程度跟他受女孩子欢迎程度是否存在一些联系啊，嘿嘿~协方差就是这样一种用来度量两个随机变量关系的统计量，我们可以仿照方差的定义：

 

来度量各个维度偏离其均值的程度，协方差可以这么来定义：

 

协方差的结果有什么意义呢？如果结果为正值，则说明两者是正相关的(从协方差可以引出“相关系数”的定义)，也就是说一个人越猥琐就越受女孩子欢迎，嘿嘿，那必须的~结果为负值就说明负相关的，越猥琐女孩子越讨厌，可能吗？如果为0，也是就是统计上说的“相互独立”。

方差是一种特殊的协方差，协方差是方差的扩展。

从协方差的定义上我们也可以看出一些显而易见的性质，如：

协方差多了就是协方差矩阵

上一节提到的猥琐和受欢迎的问题是典型二维问题，而协方差也只能处理二维问题，那维数多了自然就需要计算多个协方差，比如n维的数据集就需要计算 n! / ((n-2)!*2) 个协方差，那自然而然的我们会想到使用矩阵来组织这些数据。给出协方差矩阵的定义：

 

这个定义还是很容易理解的，我们可以举一个简单的三维的例子，假设数据集有三个维度，则协方差矩阵为

 

可见，协方差矩阵是一个对称的矩阵，而且对角线是各个维度上的方差。

Matlab协方差实战

上面涉及的内容都比较容易，协方差矩阵似乎也很简单，但实战起来就很容易让人迷茫了。必须要明确一点，协方差矩阵计算的是不同维度之间的协方差，而不是不同样本之间的。这个我将结合下面的例子说明，以下的演示将使用Matlab，为了说明计算原理，不直接调用Matlab的cov函数(蓝色部分为Matlab代码)。

首先，随机产生一个10*3维的整数矩阵作为样本集，10为样本的个数，3为样本的维数。

mysample = fix(rand(10,3)*50)

根据公式，计算协方差需要计算均值，那是按行计算均值还是按列呢，我一开始就老是困扰这个问题。前面我们也特别强调了，协方差矩阵是计算不同维度间的协方差，要时刻牢记这一点。样本矩阵的每行是一个样本，每列为一个维度，所以我们要按列计算均值。为了描述方便，我们先将三个维度的数据分别赋值：

>> dim1 = mysample(:,1);
>> dim2 = mysample(:,2);
>> dim3 = mysample(:,3);

计算dim1与dim2，dim1与dim3，dim2与dim3的协方差：

>> sum((dim1 - mean(dim1)) .* (dim2 - mean(dim2))) / (size(mysample, 1) - 1)  %得到 -147.0667
>> sum((dim1 - mean(dim1)) .* (dim3 - mean(dim3))) / (size(mysample, 1) - 1)  %得到  -82.2667
>> sum((dim2 - mean(dim2)) .* (dim3 - mean(dim3))) / (size(mysample, 1) - 1)  %得到   76.5111

搞清楚了这个后面就容易多了，协方差矩阵的对角线就是各个维度上的方差，下面我们依次计算：

>> var(dim1)  %得到 227.8778
>> var(dim2)  %得到 179.8222
>> var(dim3)  %得到 156.7111

 这样，我们就得到了计算协方差矩阵所需要的所有数据，调用Matlab自带的cov函数进行验证：

>> cov(mysample)

 把我们计算的数据对号入座，是不是一摸一样？

Update

今天突然发现，原来协方差矩阵还可以这样计算，先让样本矩阵中心化，即每一维度减去该维度的均值，使每一维度上的均值为0，然后直接用新的到的样本矩阵乘上它的转置，然后除以(N-1)即可。其实这种方法也是由前面的公式推导而来，只不过理解起来不是很直观，但在抽象的公式推导时还是很常用的！同样给出Matlab代码实现：

>> temp = mysample - repmat(mean(mysample), 10, 1);
>> result = temp' * temp ./ (size(mysample, 1) - 1)

总结

理解协方差矩阵的关键就在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间，拿到一个样本矩阵，我们最先要明确的就是一行是一个样本还是一个维度，心中明确这个整个计算过程就会顺流而下，这么一来就不会迷茫了~

下面整理了网上的代码：

1.代码：

Matlab相关系数的意义：

[cpp] view plain copy

1. Eigen::MatrixXf  correlation_matrix = corrcoef( LocM );  

对行向量求相关系数 , 与列数无关，返回 cols()*cols() 矩阵...

翻译成Eigen：
还是自己写个函数吧
//1.求协方差

[cpp] view plain copy

1. Eigen::MatrixXf CIcSearchM::cov(Eigen::MatrixXf &d1, Eigen::MatrixXf &d2)  
2. {  
3.     Eigen::MatrixXf  CovM(1,1);  
4.     assert(1 ==d1.cols() && 1 ==d2.cols() &&d1.cols()==d2.cols()  );  
5.   
6.     //求协方差  
7.     float Ex =0;float Ey=0;  
8.     for (int i=0;i< d1.rows();++i){  
9.         Ex +=d1(i);  
10.         Ey +=d2(i);  
11.     }  
12.     Ex /=d1.rows();  
13.     Ey /=d2.rows();  
14.   
15.     for (int i=0;i< d1.rows();++i){  
16.         CovM(0) += (d1(i)-Ex)*(d2(i)-Ey);  
17.     }  
18.     CovM(0) /= d1.rows() -1;  
19.     return CovM;  
20. }  

//2.写入方差矩阵

[cpp] view plain copy

1. //求矩阵的相关系数！  
2. //返回矩阵A的列向量的相关系数矩阵//对行向量求相关系数 , 与行数无关，返回 cols()*cols() 矩阵...  
3. Eigen::MatrixXf CIcSearchM::corrcoef(Eigen::MatrixXf &M)  
4. {  
5.     // C(i,j)/SQRT(C(i,i)*C(j,j)).//C is the covariation Matrix  
6.     int Row= M.rows();  
7.     int Col= M.cols();  
8.     int Order= Col;//int Order= (std::max)(Row,Col);  
9.   
10.     Eigen::MatrixXf Coef(Order,Order);  
11.     for (int i=0;i<Order;++i){  
12.         for (int j=0;j<Order;++j){  
13.             Coef(i,j)= cov((Eigen::MatrixXf)M.col(i),(Eigen::MatrixXf)M.col(j))(0);  
14.         }  
15.     }  
16.     return Coef;  
17. }  

2.优化的代码

使用Eigen计算1000维的方阵大概需要200ms的时间，相对于matlab默认开启GPU加速，时间上消耗的太多了。

参考：比较OpenBLAS、Matlab、MKL、Eigen的基础计算性能

优化的代码：

[cpp] view plain copy

1. //求矩阵的相关系数！一个原始公式的简化算法/优化算法  
2. //返回矩阵A的列向量的相关系数矩阵//对行向量求相关系数 , 与行数无关，返回 cols()*cols() 矩阵...  
3. Eigen::MatrixXf CIcSearchM::CorrcoefOpm(Eigen::MatrixXf &MI)  
4. {  
5.     Eigen::MatrixXf M =MI;  
6.     // C(i,j)/SQRT(C(i,i)*C(j,j)).//C is the covariation Matrix  
7.     //公式：  
8.     //temp = mysample - repmat(mean(mysample), 10, 1);  
9.     //result = temp' * temp ./ (size(mysample, 1) - 1)  
10.     int Row= M.rows();  
11.     int Col= M.cols();  
12.     int Order= Col;//int Order= (std::max)(Row,Col);  
13.   
14.     SYSTEMTIME sysP;   
15.     GetLocalTime( &sysP );   
16.     int MileTsp = sysP.wSecond;  
17.     int MileTP = sysP.wMilliseconds;  
18.   
19.     Eigen::MatrixXf  CovM(Order,Order);//(1,Col);  
20.     Eigen::MatrixXf  E_M(1,Col);  
21.     //减去每一个维度的均值;确定一列为一个维度。  
22.     //std::cout<< "Mat Src :"<<std::endl;m_Testor.print_EigenMat( M);  
23.     for (int i =0;i< Col;++i)  
24.     {  
25.         //求均值  
26.         E_M(i) =M.col(i).sum()/M.rows();  
27.         //std::cout<< "E_M(i)" << E_M(i)<< std::endl;  
28.         M.col(i) = M.col(i)- E_M(i);  
29.         //  
30.     }  
31.   
32.     //SYSTEMTIME sysP2;   
33.     //GetLocalTime( &sysP );   
34.     //int MileTsp2 = sysP.wSecond;  
35.     //int MileTP2 = sysP.wMilliseconds;  
36.     //int  DetaTp = MileTP2  - MileTP;  
37.     //int DetaTsp = MileTsp2 -MileTsp;  
38.     //std::cout<< "The Process time is :"<< DetaTsp<<"S"<< std::endl;  
39.     //std::cout<< "The Process time is :"<< DetaTp<<"mS"<< std::endl;  
40.   
41.     //std::cout<< "Mat E_M :"<<std::endl;m_Testor.print_EigenMat( M);  
42.     CovM = M.transpose();  
43.     //GetLocalTime( &sysP );   
44.     //MileTsp2 = sysP.wSecond;  
45.     //MileTP2 = sysP.wMilliseconds;  
46.     //DetaTp = MileTP2  - MileTP;  
47.     //DetaTsp = MileTsp2 -MileTsp;  
48.     //std::cout<< "The Process time is :"<< DetaTsp<<"S"<< std::endl;  
49.     //std::cout<< "The Process time is :"<< DetaTp<<"mS"<< std::endl;  
50.   
51.     //std::cout<< "Mat CovM :"<<std::endl;m_Testor.print_EigenMat( CovM);  
52.     CovM = CovM * M ;  
53.     //GetLocalTime( &sysP );   
54.     //MileTsp2 = sysP.wSecond;  
55.     //MileTP2 = sysP.wMilliseconds;  
56.     //DetaTp = MileTP2  - MileTP;  
57.     //DetaTsp = MileTsp2 -MileTsp;  
58.     //std::cout<< "The Process time is :"<< DetaTsp<<"S"<< std::endl;  
59.     //std::cout<< "The Process time is :"<< DetaTp<<"mS"<< std::endl;  
60.   
61.     //实现 ./ 函数 数值计算没有区别  
62.     CovM = CovM /(Order-1)/(Order-1);  
63.     //GetLocalTime( &sysP );   
64.     //MileTsp2 = sysP.wSecond;  
65.     //MileTP2 = sysP.wMilliseconds;  
66.     //DetaTp = MileTP2  - MileTP;  
67.     //DetaTsp = MileTsp2 -MileTsp;  
68.     //std::cout<< "The Process time is :"<< DetaTsp<<"S"<< std::endl;  
69.     //std::cout<< "The Process time is :"<< DetaTp<<"mS"<< std::endl;  
70.   
71.     //std::cout<< "Mat CovM :"<<std::endl;m_Testor.print_EigenMat( CovM);  
72.     //GetLocalTime( &sysP );   
73.     //MileTsp2 = sysP.wSecond;  
74.     //MileTP2 = sysP.wMilliseconds;  
75.     //DetaTp = MileTP2  - MileTP;  
76.     //DetaTsp = MileTsp2 -MileTsp;  
77.     //std::cout<< "The Process time is :"<< DetaTsp<<"S"<< std::endl;  
78.     //std::cout<< "The Process time is :"<< DetaTp<<"mS"<< std::endl;  
79.   
80.     //遍历一次  
81.     for (int i=0;i< Order;++i){  
82.         for (int j=0;j<Order;++j){  
83.             CovM(i,j) = sqrt(CovM(i,i)*CovM(j,j) );  
84.         }  
85.     }  
86.   
87.     //GetLocalTime( &sysP );   
88.     //MileTsp2 = sysP.wSecond;  
89.     //MileTP2 = sysP.wMilliseconds;  
90.     //DetaTp = MileTP2  - MileTP;  
91.     //DetaTsp = MileTsp2 -MileTsp;  
92.     //std::cout<< "The Process time is :"<< DetaTsp<<"S"<< std::endl;  
93.     //std::cout<< "The Process time is :"<< DetaTp<<"mS"<< std::endl;  
94.   
95.     //std::cout<< "Mat CovM :"<<std::endl;m_Testor.print_EigenMat( CovM);  
96.     return CovM;  
97. }