这本来是一篇时间比较久远的文章了，可是感觉既然来博客园了，也就贴出来吧。求矩阵的协方差，在很多地方都有用，我就是在用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