协方差矩阵求解算法分析

　　这本来是一篇时间比较久远的文章了，可是感觉既然来博客园了，也就贴出来吧。求矩阵的协方差，在很多地方都有用，我就是在用Matlab做数字图像处理时用到的这个。为了理解，我看了一下午的书，什么线性代数，什么概率论都被我翻出来了，把思路贴一下吧：

一、定义

设n维随机变量(X1,X2,...,Xn)的协方差
         c(i,j) = cov(X(i),X(j)) = E{[X(i)-E(X(i))][X(j)-E(X(j))]}      i,j=1,2,...,n   （E是期望，即平均值）
都存在，则称矩阵
      c(1,1)  c(1,2)  ...  c(1,n)
      c(2,1)  c(2,2)  ...  c(2,n)
C =      .       .           .
            .       .           .
      c(n,1)  c(n,2)  ...  c(n,n)
为n维随机变量(X1,X2,...,Xn)的协方差矩阵，由于c(i,j)=c(j,i)，因而上述矩阵是一个对称矩阵。

 

二、矩阵的协方差

        X为矩阵，则它的每一列被视为一个向量，如果是一个9 X 3的矩阵，即是一个三维的向量（列1，列2，列3）。协方差矩阵的元素为c(i,j) = E{[列i-E(列i)]*[列j-E(列j)]}。
        在Matlab环境下设计算法有：
[n,p] = size(X); 
方法一：c = cov(X,1);（自带的- -）
方法二：for   i=1:p
                 for   j=1:p
                    c(i,j)=mean((X(:,i)-mean(X(:,i))).*(X(:,j)-mean(X(:,j))));
                 end
              end   
方法三：Y = X-ones(n,1)*mean(X); 
              c = Y'*Y/n; 
方法四：c = X'*X/n-mean(X)'*mean(X);
        这里需要说明的是，matlab中计算协方差是首先加了n项，然后除以n或n-1，cov(x)呢是除以n-1，cov(x,1)呢是除以n。cov(x)是无偏估计,cov(x,1)是最大似然估计,cov(x,1)=cov(x)*(n-1)/n。我在这里计算的都是最大似然估计。

 

三、测试与改进
        我定义了一个矩阵X = [1 1 2;-2 3 1;4 0 3]，然后用X = repmat(X,3000000,1)竖向复制了3,000,000份，组成一个9000000×3的测试矩阵。
        经过测试结果都一致：
c =
    6.8889   -2.7778    2.0000
   -2.7778    1.5556   -1.0000
    2.0000   -1.0000    0.6667

        可改进算法，将mean(X)，即平均值那一项，提取出来：
[n,p] = size(X); 
avg = mean(X);
方法二：for   i=1:p
                 for   j=1:p
                     c(i,j)=mean((X(:,i)-avg(i)).*(X(:,j)-avg(j)));
                 end
              end   
方法三：Y = X-ones(n,1)*avg; 
              c = Y'*Y/n; 
方法四：c = X'*X/n-avg'*avg;
        用时如下：
方法一用时1.516000 seconds.
方法二原用时8.891000 seconds，改进后为6.125000 seconds.
方法三原用时1.703000 seconds，改进后为1.547000 seconds.
方法四原用时0.891000 seconds，改进后为0.547000 seconds.
        以上是在我的机器上测试，配置CPU为AMD 4400+ 2.2GHz X2，2GB内存。可以看出，向量化循环可以大大的减少计算时间啊！

 

四、小细节

      好吧，就是这样了，其它的小改动也有，但是算法基本就是三个，这三个算法的细节如果做调整，如方法三的Y，如果直接X=X-ones(n,1)*avg的话，可以节省空间成本，时间成本也会少0.1秒不到的样子，但是会破坏原始数据；如果把ones(n,1)*avg改成repmat(avg,n,1)的话，时间还会增加0.2秒以上；如果想求无偏估计的话，把方法三、四只需要把最后除的n改成n-1即可，方法二则要改成c(i,j)=sum((X(:,i)-avg(i)).*(X(:,j)-avg(j)))/(n-1)。

转载请注明原址：http://www.cnblogs.com/lekko/archive/2012/07/20/2600966.html