Ketidakstabilan numerik menghitung matriks kovarians terbalik

Saya memiliki 65 sampel data 21 dimensi (disisipkan di sini ) dan saya membuat matriks kovarians. Ketika dihitung dalam C ++ saya mendapatkan matriks kovarians disisipkan di sini . Dan ketika dihitung dalam matlab dari data (seperti yang ditunjukkan di bawah ini) saya mendapatkan matriks kovarians yang disisipkan di sini

Kode Matlab untuk menghitung cov dari data:

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

Seperti yang Anda lihat perbedaan dalam matriks kovarian adalah menit (~ e-07), yang mungkin disebabkan oleh masalah numerik dalam kompiler menggunakan aritmatika floating point.

Namun, ketika saya menghitung matriks kovarians semu-terbalik dari matriks kovarians yang diproduksi oleh matlab dan yang dihasilkan oleh kode C ++ saya, saya mendapatkan hasil yang sangat berbeda. Saya menghitungnya dengan cara yang sama yaitu:

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);

Perbedaannya sangat besar sehingga ketika saya menghitung jarak mahalanobis dari sampel (disisipkan di sini ) ke distribusi 65 sampel dengan:

$(65/64^2) \times ((sample-mean)\times {\sum}^{-1} \times (sample-mean)')$

menggunakan matriks kovarians terbalik yang berbeda ( ) Saya mendapatkan hasil yang sangat berbeda yaitu:${\sum}^{-1}$

 (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

Apakah normal jika perbedaan kecil (e-7) dalam matriks kovarians memiliki efek pada perhitungan matriks pseudo-invers? Dan jika demikian, apa yang bisa saya lakukan untuk mengurangi efek ini?

Gagal ini, apakah ada metrik jarak lain yang dapat saya gunakan yang tidak melibatkan kovarians terbalik? Saya menggunakan jarak Mahalanobis seperti yang kita tahu untuk n sampel itu mengikuti distribusi beta, yang saya gunakan untuk pengujian hipotesis

Banyak terima kasih sebelumnya

EDIT: Menambahkan C ++ kode untuk menghitung kovarians matriks di bawah ini: The vector<vector<double> >mewakili koleksi baris dari file disisipkan.

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; 
        }
    }

Matriks yang Anda cari untuk dibalik bukanlah matriks kovarian yang "valid" karena mereka tidak pasti positif; secara numerik mereka bahkan memiliki beberapa nilai eigen yang negatif (tetapi mendekati nol). Ini kemungkinan besar disebabkan oleh nol mesin, misalnya nilai eigen terakhir dari matriks "matlab_covariance" Anda adalah -0.000000016313723. Untuk mengoreksi ke definitif positif, Anda dapat melakukan dua hal:

1. Tambahkan saja konstanta yang sangat kecil di sepanjang diagonal; bentuk koreksi Tikhonov benar-benar ( dengan ).$Χ+\lambda I$$\lambda \rightarrow 0$
2. (Berdasarkan apa yang diusulkan whuber) Gunakan SVD, atur nilai eigen "bermasalah" menjadi nilai kecil tetap (bukan nol), merekonstruksi matriks kovarians Anda dan kemudian membalikkannya. Jelas jika Anda menetapkan beberapa nilai eigen tersebut menjadi nol, Anda akan berakhir dengan matriks non-negatif (atau semi-positif), yang tidak akan dapat dibalik lagi.

Matriks non-negatif tidak memiliki invers tetapi ia memiliki invers pseudo (semua matriks dengan entri real atau kompleks memiliki invers semu), namun invers pseudo-invers Moore-Penrose lebih mahal secara komputasi daripada invers yang benar dan jika invers yang ada itu sama dengan invers semu. Jadi pergi saja untuk kebalikannya :)

Kedua metode praktis mencoba untuk menangani nilai eigen yang mengevaluasi nol (atau di bawah nol). Metode pertama agak bergelombang tetapi mungkin jauh lebih cepat untuk diterapkan. Untuk sesuatu yang sedikit lebih stabil, Anda mungkin ingin menghitung SVD dan kemudian mengatur sama dengan absolut dari nilai eigen terkecil (sehingga Anda mendapatkan non-negatif) ditambah sesuatu yang sangat kecil (sehingga Anda mendapatkan positif). Hanya berhati-hatilah untuk tidak memaksakan kepositifan ke matriks yang jelas negatif (atau sudah positif). Kedua metode akan mengubah nomor pengkondisian matriks Anda.$\lambda$

Dalam istilah statistik apa yang Anda lakukan dengan menambahkan di diagonal matriks kovarian Anda, Anda menambahkan noise ke pengukuran Anda. (Karena diagonal dari matriks kovarians adalah varian dari setiap titik dan dengan menambahkan sesuatu ke nilai-nilai itu, Anda hanya mengatakan "varians pada titik yang saya baca sebenarnya sedikit lebih besar daripada yang saya pikir awalnya".)$\lambda$

Tes cepat untuk kepastian positif dari sebuah matriks adalah keberadaan (atau tidak) dekomposisi Cholesky itu.

Juga sebagai catatan komputasi:

1. Jangan gunakan pelampung, gunakan ganda untuk menyimpan matriks kovarians Anda.
2. Gunakan pustaka aljabar linier numerik dalam C ++ (seperti Eigen atau Armadillo) untuk mendapatkan inversi dari matriks, produk matriks, dll. Lebih cepat, lebih aman, dan lebih ringkas.

EDIT: Mengingat Anda memiliki dekomposisi Cholesky dari matriks Anda sehingga (Anda harus melakukan itu untuk memeriksa Anda memiliki matriks Pos.Def.) Anda harus dapat segera menyelesaikan sistem . Anda baru saja menyelesaikan Ly = b untuk y dengan substitusi maju, dan kemudian L ^ Tx = y untuk x dengan substitusi kembali. (Dalam eigen cukup gunakan .solve (x) metode objek Cholesky Anda) Terima kasih kepada bnaul dan Zen untuk menunjukkan bahwa saya sangat fokus untuk mendapatkan be Pos.Def. bahwa saya lupa mengapa kami peduli tentang itu di tempat pertama :)$K$$LL^T$$Kx=b$$K$

Jawaban dan komentar yang diposting semua memberikan poin bagus tentang bahaya membalikkan matriks yang hampir tunggal. Namun, sejauh yang saya tahu, tidak ada yang menyebutkan bahwa menghitung jarak Mahalanobis sebenarnya tidak memerlukan pembalikan kovarians sampel. Lihat pertanyaan StackOverflow ini untuk deskripsi bagaimana melakukannya menggunakan dekomposisi .$LU$

Prinsipnya sama dengan menyelesaikan sistem linier: ketika mencoba memecahkan untuk sehingga , ada jauh lebih efisien dan metode yang stabil secara numerik daripada mengambil .$x$$Ax=b$$x=A^{-1}b$

Sunting: mungkin tidak perlu dikatakan, tetapi metode ini menghasilkan nilai jarak yang tepat, sedangkan menambahkan ke dan membalikkan menghasilkan hanya perkiraan.$\lambda I$$S$

(Bertahun-tahun kemudian) contoh kecil: dengan peringkat-kekurangan, eigen dari akan 0 dalam mesin presisi - dan sekitar setengah dari "nol" mungkin :$A$$r < n, \ n - r$$A^T A$$< 0$

#!/usr/bin/env python2
""" many eigenvalues of A'A are tiny but < 0 """
# e.g. A 1 x 10: [-1.4e-15 -6.3e-17 -4e-17 -2.7e-19 -8.8e-21  1e-18 1.5e-17 5.3e-17 1.4e-15  7.7]

from __future__ import division
import numpy as np
from numpy.linalg import eigvalsh  # -> lapack_lite
# from scipy.linalg import eigvalsh  # ~ same
from scipy import __version__

np.set_printoptions( threshold=20, edgeitems=10, linewidth=140,
        formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
print "versions: numpy %s  scipy %s \n" % (
        np.__version__, __version__  )

np.random.seed( 3 )

rank = 1
n = 10
A = np.random.normal( size=(rank, n) )
print "A: \n", A
AA = A.T.dot(A)
evals = eigvalsh( AA )
print "eigenvalues of A'A:", evals