C++ PCA - 计算协方差矩阵

Question

我正在尝试用 C++ 计算 vector 的协方差(矩阵)...

我已经执行了以下操作:

std::vector<std::vector<double> > data = { {2.5, 2.4}, {0.5, 0.7} };

然后我计算并减去平均值，得到以下结果:

data = { {0.05, -0.05}, {-0.1, 0.1} }

据我所知，下一步是转置矩阵，将原点相乘，求和，最后除以维度 X - 1..

我写了以下内容:

void cover(std::vector<std::vector<double> > &d)
{
    double cov = 0.0; 

    for(unsigned i=0; (i < d.size()); i++)
    {
        for(unsigned j=0; (j < d[i].size()); j++)
        {
            cov += d[i][j] * d[j][i] / (d[i].size() - 1);
            std::cout << cov << " ";
        }
        std::cout << std::endl;
    }
}

其中 d 是从每个点中减去平均值后的 vector 这给了我结果:

0.0025, 0.0075 
0.0125, 0.0225

与matlab相比:

2.0000    1.7000
1.7000    1.4450

有人知道我哪里出错了吗？

谢谢

Answer 1

这个声明:

As far as I'm aware, the next step is to transpose the matrix, and multiply the origin together, take the sum and finally divide by the dimensions X - 1..

这个实现:

cov += d[i][j] * d[j][i] / (d[i].size() - 1);

不要说同样的话。基于定义here :

void outer_product(vector<double> row, vector<double> col, vector<vector<double>>& dst) {
    for(unsigned i = 0; i < row.size(); i++) {
        for(unsigned j = 0; j < col.size(); i++) {
            dst[i][j] = row[i] * col[j];
        }
    }
}

//computes row[i] - val for all i;
void subtract(vector<double> row, double val, vector<double>& dst) {
    for(unsigned i = 0; i < row.size(); i++) {
        dst[i] = row[i] - val;
    }
}

//computes m[i][j] + m2[i][j]
void add(vector<vector<double>> m, vector<vector<double>> m2, vector<vector<double>>& dst) {
    for(unsigned i = 0; i < m.size(); i++) {
        for(unsigned j = 0; j < m[i].size(); j++) { 
            dst[i][j] = m[i][j] + m2[i][j];
        }
    }
}

double mean(std::vector<double> &data) {
    double mean = 0.0;

    for(unsigned i=0; (i < data.size());i++) {
        mean += data[i];
    }

    mean /= data.size();
    return mean;
}

void scale(vector<vector<double>> & d, double alpha) {
    for(unsigned i = 0; i < d.size(); i++) {
        for(unsigned j = 0; j < d[i].size(); j++) {
            d[i][j] *= alpha;
        }
    }
}

因此，根据这些定义，我们可以计算协方差矩阵的值。

void compute_covariance_matrix(vector<vector<double>> & d, vector<vector<double>> & dst) {
    for(unsigned i = 0; i < d.size(); i++) {
        double y_bar = mean(d[i]);
        vector<double> d_d_bar(d[i].size());
        subtract(d[i], y_bar, d_d_bar);
        vector<vector<double>> t(d.size());
        outer_product(d_d_bar, d_d_bar, t);
        add(dst, t, dst);
    } 
    scale(dst, 1/(d.size() - 1));
}