Eigen: существует ли встроенный способ вычисления выборочной ковариации

Я использую библиотеку Eigen в С++: в настоящее время я вычисляю матрицу ковариации следующим образом:

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

Есть ли встроенный/более оптимизированный способ сделать это с помощью библиотеки Eigen? Например, если я храню свои данные в MatrixXd, где каждая строка является наблюдением, а каждый столбец - функцией?

Спасибо

Ответы

Ответ 1

Когда каждая строка является наблюдением, вы можете использовать матричную формулу для матрицы ковариационной выборки, как показано на wikipedia (http://en.wikipedia.org/wiki/Sample_mean_and_sample_covariance#Sample_covariance)

Sample covariance, source: wikipedia article linked above.

Это довольно легко написать в терминах собственных умножений матрицы и т.д. Является ли это более показательным для меня, не очевидно, я подозреваю, что оптимизатор должен будет сделать действительно хорошую работу (обязательно используйте хотя бы - O2). Возможно, стоит попробовать и профилировать его.

Ответ 2

Использование собственных выражений будет использовать алгоритмы с оптимизацией SIMD и кешей, поэтому да, это определенно будет быстрее и в любом случае гораздо проще написать:

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

Кроме того, если предположить, что "данные" являются typedef для double [21], то вы можете использовать функцию Map < > для просмотра вашего std::vector как объекта Eigen:

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