計算逆協方差矩陣的數值不穩定性

我有 65 個 21 維數據樣本（粘貼在此處），我正在從中構建協方差矩陣。當用 C++ 計算時，我得到了粘貼在這裡的協方差矩陣。當在 matlab 中從數據中計算時（如下所示），我得到了粘貼在這裡的協方差矩陣

用於從數據計算 cov 的 Matlab 代碼：

data = csvread('path/to/data');
matlab_cov = cov(data);

如您所見，協方差矩陣的差異很小（~e-07），這可能是由於使用浮點運算的編譯器中的數值問題。

但是，當我從 matlab 生成的協方差矩陣和我的 C++ 代碼生成的協方差矩陣計算偽逆協方差矩陣時，我得到了截然不同的結果。我以相同的方式計算它們，即：

data = csvread('path/to/data');
matlab_cov = cov(data);
my_cov = csvread('path/to/cov_file');
matlab_inv = pinv(matlab_cov);
my_inv = pinv(my_cov);

差異是如此之大，以至於當我計算從樣本（粘貼在此處）到 65 個樣本的分佈的馬氏距離時：

使用不同的逆協方差矩陣（) 我得到了截然不同的結果，即：

(65/(64^2))*((sample-sample_mean)*my_inv*(sample-sample_mean)')
ans =

  1.0167e+05

(65/(64^2))*((sample-sample_mean)*matlab_inv*(sample-sample_mean)')
ans =

 109.9612

協方差矩陣中的小（e-7）差異對偽逆矩陣的計算產生這種影響是否正常？如果是這樣，我能做些什麼來減輕這種影響？

如果做不到這一點，我可以使用其他不涉及逆協方差的距離度量嗎？我使用馬氏距離，因為我們知道 n 個樣本遵循 beta 分佈，我用它來進行假設檢驗

提前謝謝了

編輯：在下面添加用於計算協方差矩陣的 C++ 代碼： 表示vector<vector<double> >粘貼文件中的行集合。

Mat covariance_matrix = Mat(21, 21, CV_32FC1, cv::Scalar(0));
   for(int j = 0; j < 21; j++){
       for(int k = 0; k < 21; k++){
           for(std::vector<vector<double> >::iterator it = data.begin(); it!= data.end(); it++){
               covariance_matrix.at<float>(j,k) += (it->at(j) - mean.at(j)) * (it->at(k) - mean[k]);
           }
           covariance_matrix.at<float>(j,k) /= 64; 
       }
   }

您要反轉的矩陣不是“有效的”協方差矩陣，因為它們不是正定的；從數值上看，它們甚至有一些負的特徵值（但接近於零）。這很可能是由於機器零，例如“matlab_covariance”矩陣的最後一個特徵值是-0.000000016313723。要更正為肯定，您可以做兩件事：

1. 只需沿對角線添加一些非常小的常數；一種形式的 Tikhonov 校正真的（和).
2. （基於 whuber 的建議）使用 SVD，將“有問題的”特徵值設置為固定的小值（非零），重建協方差矩陣，然後將其反轉。顯然，如果您將其中一些特徵值設置為零，您最終將得到一個非負（或半正）矩陣，該矩陣仍然不可逆。

非負矩陣沒有逆矩陣，但它確實有一個偽逆矩陣（所有具有實數或複數元素的矩陣都有一個偽逆矩陣），但是 Moore-Penrose 偽逆矩陣在計算上比真正的逆矩陣更昂貴，並且如果逆存在它等於偽逆。所以就反過來吧:)

這兩種方法實際上都試圖處理評估為零（或低於零）的特徵值。第一種方法有點手動，但實施起來可能要快得多。對於稍微更穩定的東西，您可能需要計算 SVD，然後設置等於最小特徵值的絕對值（所以你得到非負數）加上非常小的東西（所以你得到正數）。請注意不要對明顯為負（或已經為正）的矩陣強制執行正數。這兩種方法都會改變矩陣的條件數。

在統計術語中，您通過添加在協方差矩陣的對角線上，您可以在測量中添加噪聲。（因為協方差矩陣的對角線是每個點的方差，並且通過向這些值添加一些內容，您只需說“我所讀取的點的方差實際上比我最初認為的要大一些”。）

一個矩陣正定性的快速測試是它的 Cholesky 分解的存在（或不存在）。

也作為計算說明：

1. 不要使用浮點數，使用雙精度數來存儲協方差矩陣。
2. 使用 C++ 中的數值線性代數庫（如 Eigen 或 Armadillo）來獲得矩陣的逆、矩陣乘積等。它更快、更安全、更簡潔。

編輯：假設您對矩陣進行了 Cholesky 分解這樣（您必須這樣做以檢查您是否擁有 Pos.Def. 矩陣）您應該能夠立即解決系統問題. 您只需通過前向替換解決 Ly = b 為 y，然後通過反向替換為 x 解決 L^Tx = y。（在 eigen 中，只需使用 Cholesky 對象的 .solve(x) 方法）感謝 bnaul 和 Zen 指出我非常專注於獲取成為Pos.Def。我忘記了為什麼我們首先關心這個:)

引用自：https://stats.stackexchange.com/questions/50947