共變數矩陣求解演算法分析

　　這本來是一篇時間比較久遠的文章了，可是感覺既然來部落格園了，也就貼出來吧。求矩陣的共變數，在很多地方都有用，我就是在用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