天天看點

numpy和MATLAB中協方差矩陣的計算

前言:最近重新看了下PCA算法的實作,在裡面用到了協方差矩陣的計算。為了今後友善回顧,特此記錄。

這裡隻針對協方差矩陣的計算做說明,并與numpy和MATLAB中的協方差矩陣計算做對比。關于協方差矩陣的定義這裡不再贅述。

混淆點

關于矩陣中資料是按照列排列和按照行排列的出來的協方差矩陣是不同的。

如果資料按照列排列,例如一個3*2的矩陣,那麼該矩陣有2個資料(或者叫觀測值),3個變量。

如果資料按照行排列,同樣以3*2的矩陣為例,那麼該矩陣有3個資料(或者叫觀察值),2個變量。

不過不管矩陣中資料是按照行排列還是列排列,最後的協方差矩陣的次元隻與變量有關。

計算

下面将分别解釋資料按行排列和列排列的協方差矩陣的計算。

一、資料按照行排列

矩陣A包含2個觀察值,每個觀察值有3個變量。

numpy和MATLAB中協方差矩陣的計算

①計算變量的均值(也就是每列的均值)

numpy和MATLAB中協方差矩陣的計算

 ②變量減去均值

numpy和MATLAB中協方差矩陣的計算

③計算協方差矩陣(m代表資料的個數)

numpy和MATLAB中協方差矩陣的計算

 ④代碼驗證

numpy

In [39]: A = np.array([[1,5,3],[6,7,2]])

In [40]: Amean = np.mean(A,axis=0)

In [41]: U = A - Amean

In [42]: U.T @ U /(2-1)
Out[42]:
array([[12.5,  5. , -2.5],
       [ 5. ,  2. , -1. ],
       [-2.5, -1. ,  0.5]])

In [43]: np.cov(A,rowvar=0)#numpy中預設是行是變量,列是資料。所有将rowvar置0,代表行為資料
Out[43]:
array([[12.5,  5. , -2.5],
       [ 5. ,  2. , -1. ],
       [-2.5, -1. ,  0.5]])
           

MATLAB

>> A = [1 5 3;6 7 2]

A =

     1     5     3
     6     7     2

>> cov(A)

ans =

   12.5000    5.0000   -2.5000
    5.0000    2.0000   -1.0000
   -2.5000   -1.0000    0.5000
           

二、資料按照列排列

矩陣A包含3個觀察值,每個觀察值有2個變量。

numpy和MATLAB中協方差矩陣的計算

①計算變量的均值(也就是每行的均值)

numpy和MATLAB中協方差矩陣的計算

 ②變量減去均值

numpy和MATLAB中協方差矩陣的計算

 ③計算協方差矩陣(m代表資料的個數)

numpy和MATLAB中協方差矩陣的計算

 ④代碼驗證

In [53]: A = np.array([[1,5,3],[6,7,2]])

In [54]: Amean = np.mean(A,axis=1,keepdims=True)

In [55]: U = A - Amean

In [56]: U @ U.T / (3 - 1)
Out[56]:
array([[4., 1.],
       [1., 7.]])

In [57]: np.cov(A)#numpy預設行是變量
Out[57]:
array([[4., 1.],
       [1., 7.]])
           

總結

在計算矩陣的協方差矩陣時關鍵點是要弄清楚矩陣是一行對應一組資料還是一列對應一組資料,計算出來的協方差矩陣都是n維的方陣(n代表變量的個數)。在計算時,第一步都是用原始矩陣減去變量的均值得到變量和為0的矩陣U;第二步根據資料的排列方式進行計算,若一行代表一個資料則選用公式U.T*U/(m-1),若一列代表一個資料則U*U/(m-1)。