用c++实现pca算法

PCA（Principal Component Analysis）是一种常用的数据降维方法，可以将高维数据映射到低维空间中，同时保留数据的主要特征。下面是用C++实现PCA算法的示例代码：

#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

// 计算矩阵的转置
vector<vector<double>> transpose(vector<vector<double>> matrix) {
    int rows = matrix.size();
    int cols = matrix[0].size();
    vector<vector<double>> res(cols, vector<double>(rows));
    for(int i = 0; i < cols; i++) {
        for(int j = 0; j < rows; j++) {
            res[i][j] = matrix[j][i];
        }
    }
    return res;
}

// 计算矩阵的乘积
vector<vector<double>> dot(vector<vector<double>> A, vector<vector<double>> B) {
    int rowsA = A.size();
    int colsA = A[0].size();
    int rowsB = B.size();
    int colsB = B[0].size();
    vector<vector<double>> res(rowsA, vector<double>(colsB));
    for(int i = 0; i < rowsA; i++) {
        for(int j = 0; j < colsB; j++) {
            double sum = 0.0;
            for(int k = 0; k < colsA; k++) {
                sum += A[i][k] * B[k][j];
            }
            res[i][j] = sum;
        }
    }
    return res;
}

// 计算矩阵的特征值和特征向量
void eigen(vector<vector<double>> matrix, vector<double>& eigenvalues, vector<vector<double>>& eigenvectors) {
    int n = matrix.size();
    vector<vector<double>> A = matrix;
    vector<vector<double>> V(n, vector<double>(n));
    for(int i = 0; i < n; i++) {
        V[i][i] = 1.0;
    }
    int max_iterations = 1000; // 最大迭代次数
    double eps = 1e-6; // 精度控制
    int iter = 0;
    while(iter < max_iterations) {
        double sm = 0.0; // 计算非对角线元素的平方和
        for(int i = 0; i < n - 1; i++) {
            for(int j = i + 1; j < n; j++) {
                sm += A[i][j] * A[i][j];
            }
        }
        if(sm < eps) {
            break; // 达到精度要求
        }
        for(int p = 0; p < n - 1; p++) {
            for(int q = p + 1; q < n; q++) {
                if(A[p][q] != 0) {
                    double app = A[p][p];
                    double aqq = A[q][q];
                    double apq = A[p][q];
                    double phi = 0.5 * atan2(2 * apq, aqq - app); // 计算旋转角度
                    double c = cos(phi);
                    double s = sin(phi);
                    double app1 = c * c * app - 2 * s * c * apq + s * s * aqq;
                    double aqq1 = s * s * app + 2 * s * c * apq + c * c * aqq;
                    if(abs(aqq1 - aqq) > eps && abs(app1 - app) > eps) {
                        A[p][p] = app1;
                        A[q][q] = aqq1;
                        A[p][q] = 0.0;
                        A[q][p] = 0.0;
                        for(int i = 0; i < p; i++) {
                            double api = A[i][p];
                            double aqi = A[i][q];
                            A[i][p] = c * api - s * aqi;
                            A[i][q] = s * api + c * aqi;
                        }
                        for(int i = p + 1; i < q; i++) {
                            double api = A[p][i];
                            double aqi = A[i][q];
                            A[p][i] = c * api - s * aqi;
                            A[i][q] = s * api + c * aqi;
                        }
                        for(int i = q + 1; i < n; i++) {
                            double api = A[p][i];
                            double aqi = A[q][i];
                            A[p][i] = c * api - s * aqi;
                            A[q][i] = s * api + c * aqi;
                        }
                        for(int i = 0; i < n; i++) {
                            double vpi = V[i][p];
                            double vqi = V[i][q];
                            V[i][p] = c * vpi - s * vqi;
                            V[i][q] = s * vpi + c * vqi;
                        }
                    }
                }
            }
        }
        iter++;
    }
    for(int i = 0; i < n; i++) {
        eigenvalues.push_back(A[i][i]);
        eigenvectors.push_back(V[i]);
    }
}

// 计算PCA
vector<vector<double>> pca(vector<vector<double>> data, int k) {
    int n = data.size();
    int m = data[0].size();
    vector<vector<double>> X = data;
    // 中心化处理
    vector<double> mean(m);
    for(int j = 0; j < m; j++) {
        double sum = 0.0;
        for(int i = 0; i < n; i++) {
            sum += X[i][j];
        }
        mean[j] = sum / n;
        for(int i = 0; i < n; i++) {
            X[i][j] -= mean[j];
        }
    }
    // 计算协方差矩阵
    vector<vector<double>> Cov = dot(transpose(X), X);
    Cov = Cov / (n - 1);
    // 计算特征值和特征向量
    vector<double> eigenvalues;
    vector<vector<double>> eigenvectors;
    eigen(Cov, eigenvalues, eigenvectors);
    // 取前k个特征向量
    vector<vector<double>> W(k, vector<double>(m));
    for(int j = 0; j < k; j++) {
        for(int i = 0; i < m; i++) {
            W[j][i] = eigenvectors[j][i];
        }
    }
    // 将数据映射到低维空间中
    vector<vector<double>> Y = dot(X, transpose(W));
    return Y;
}

int main() {
    vector<vector<double>> data = {{2.5, 2.4}, {0.5, 0.7}, {2.2, 2.9}, {1.9, 2.2}, {3.1, 3.0}, {2.3, 2.7}, {2.0, 1.6}, {1.0, 1.1}, {1.5, 1.6}, {1.1, 0.9}};
    vector<vector<double>> Y = pca(data, 1);
    for(int i = 0; i < Y.size(); i++) {
        cout << Y[i][0] << endl;
    }
    return 0;
}

上面的代码中，`transpose`函数用于计算矩阵的转置，`dot`函数用于计算矩阵的乘积，`eigen`函数用于计算矩阵的特征值和特征向量，`pca`函数用于计算PCA。

注：上述代码中的运算符“/”表示矩阵的逐元素除法。