Eigen: Czy jest wbudowana sposób obliczyć przykładowy kowariancji

Question

używam biblioteki EIGEN w C++: Jestem obecnie obliczania macierzy kowariancji się następująco:

Eigen::MatrixXd covariance_matrix = Eigen::MatrixXd::Constant(21, 21, 0); 
data mean = calc_mean(all_data) 
for(int j = 0; j < 21; j++){ 
    for(int k = 0; k < 21; k++){ 
     for(std::vector<data>::iterator it = all_data.begin(); it!= all_data.end(); it++){ 
      covariance_matrix(j,k) += ((*it)[j] - mean[j]) * ((*it)[k] - mean[k]); 
     } 
     covariance_matrix(j,k) /= all_data.size() - 1; 
    } 
} 

Czy jest wbudowana/bardziej zoptymalizowany sposób zrobić to z biblioteką Eigen? Na przykład, jeśli przechowuję moje dane w MatrixXd, gdzie każdy wiersz jest obserwacją, a każda kolumna jest funkcją?

Dzięki

Answer 1

Gdy każdy rząd jest obserwacja, można użyć preparatu matrycowego o próbki macierzy kowariancji jak pokazano na Wikipedia (http://en.wikipedia.org/wiki/Sample_mean_and_sample_covariance#Sample_covariance)

.

Jest to dość łatwe do napisania pod względem mnożenia macierzy Eigen itd. Bez względu na to, czy będzie to bardziej wydajne, nie jest dla mnie oczywiste, podejrzewam, że optymalizator musiałby wykonać naprawdę dobrą robotę (pamiętaj, aby użyć co najmniej -O2). Warto spróbować i profilować go.

Answer 2

Korzystanie z wyrażeń EIGEN wykorzysta SIMD oraz algorytmy cache zoptymalizowany, więc tak powinno zdecydowanie szybciej, aw każdym razie znacznie prostsze napisać:

MatrixXd centered = mat.rowwise() - mat.colwise().mean(); 
MatrixXd cov = (centered.adjoint() * centered)/double(mat.rows() - 1); 

Ponadto, zakładając, że „dane” jest typedef dla podwójny [21], a następnie można użyć funkcji> Mapa < aby wyświetlić std :: vector jako obiekt Eigen:

Map<Matrix<double,Dynamic,21,RowMajor> > mat(&(all_data[0][0], all_data.size(), 21);