c语言用鸢尾花数据集实现pca

PCA（Principal Component Analysis）是一种常用的数据降维算法，可以将高维数据降到低维空间，以便于可视化和处理。

使用鸢尾花数据集实现PCA的步骤如下：

1. 导入头文件和数据集

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 150  // 样本数
#define D 4    // 特征数

double data[N][D] = {0};  // 数据集

2. 读取数据集

int read_data(const char* filename) {
    FILE* fp = fopen(filename, "r");
    if (fp == NULL) {
        printf("Failed to open file: %s\n", filename);
        return -1;
    }

    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < D; ++j) {
            fscanf(fp, "%lf", &data[i][j]);
        }
    }

    fclose(fp);
    return 0;
}

3. 中心化数据集

void center_data() {
    double mean[D] = {0};  // 每个特征的均值
    for (int j = 0; j < D; ++j) {
        for (int i = 0; i < N; ++i) {
            mean[j] += data[i][j];
        }
        mean[j] /= N;
    }

    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < D; ++j) {
            data[i][j] -= mean[j];
        }
    }
}

4. 计算协方差矩阵

void cov(double* covmat) {
    for (int i = 0; i < D; ++i) {
        for (int j = 0; j < D; ++j) {
            double sum = 0;
            for (int k = 0; k < N; ++k) {
                sum += data[k][i] * data[k][j];
            }
            covmat[i*D+j] = sum / (N-1);
        }
    }
}

5. 计算特征值和特征向量

void eig(double* covmat, double* eigval, double* eigvec) {
    // 使用Jacobi方法求解特征值和特征向量
    // 可以参考 https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
    // 这里为了简单起见，直接使用了Eigen库
    // 如果需要自己实现Jacobi方法，可以参考 https://github.com/yixuan/RobustPCA/blob/master/src/jacobi.cpp
    // 或者 https://github.com/SolbiatiAlessandro/CppNumericalSolvers/blob/main/src/eigen/jacobi_eigensolver.cpp
    Eigen::Map<Eigen::MatrixXd> cov(covmat, D, D);
    Eigen::Map<Eigen::VectorXd> val(eigval, D);
    Eigen::Map<Eigen::MatrixXd> vec(eigvec, D, D);
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
    val = eig.eigenvalues();
    vec = eig.eigenvectors();
}

6. 取前k个主成分

void pca(int k, double* pc) {
    center_data();

    double covmat[D*D] = {0};
    cov(covmat);

    double eigval[D] = {0};
    double eigvec[D*D] = {0};
    eig(covmat, eigval, eigvec);

    for (int i = 0; i < k; ++i) {
        for (int j = 0; j < D; ++j) {
            pc[i*D+j] = eigvec(j, D-1-i);
        }
    }
}

完整代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <Eigen/Eigen>

#define N 150  // 样本数
#define D 4    // 特征数

double data[N][D] = {0};  // 数据集

int read_data(const char* filename) {
    FILE* fp = fopen(filename, "r");
    if (fp == NULL) {
        printf("Failed to open file: %s\n", filename);
        return -1;
    }

    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < D; ++j) {
            fscanf(fp, "%lf", &data[i][j]);
        }
    }

    fclose(fp);
    return 0;
}

void center_data() {
    double mean[D] = {0};  // 每个特征的均值
    for (int j = 0; j < D; ++j) {
        for (int i = 0; i < N; ++i) {
            mean[j] += data[i][j];
        }
        mean[j] /= N;
    }

    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < D; ++j) {
            data[i][j] -= mean[j];
        }
    }
}

void cov(double* covmat) {
    for (int i = 0; i < D; ++i) {
        for (int j = 0; j < D; ++j) {
            double sum = 0;
            for (int k = 0; k < N; ++k) {
                sum += data[k][i] * data[k][j];
            }
            covmat[i*D+j] = sum / (N-1);
        }
    }
}

void eig(double* covmat, double* eigval, double* eigvec) {
    // 使用Jacobi方法求解特征值和特征向量
    // 可以参考 https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
    // 这里为了简单起见，直接使用了Eigen库
    // 如果需要自己实现Jacobi方法，可以参考 https://github.com/yixuan/RobustPCA/blob/master/src/jacobi.cpp
    // 或者 https://github.com/SolbiatiAlessandro/CppNumericalSolvers/blob/main/src/eigen/jacobi_eigensolver.cpp
    Eigen::Map<Eigen::MatrixXd> cov(covmat, D, D);
    Eigen::Map<Eigen::VectorXd> val(eigval, D);
    Eigen::Map<Eigen::MatrixXd> vec(eigvec, D, D);
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
    val = eig.eigenvalues();
    vec = eig.eigenvectors();
}

void pca(int k, double* pc) {
    center_data();

    double covmat[D*D] = {0};
    cov(covmat);

    double eigval[D] = {0};
    double eigvec[D*D] = {0};
    eig(covmat, eigval, eigvec);

    for (int i = 0; i < k; ++i) {
        for (int j = 0; j < D; ++j) {
            pc[i*D+j] = eigvec(j, D-1-i);
        }
    }
}

int main() {
    if (read_data("iris.data") < 0) {
        return -1;
    }

    double pc[2*D] = {0};
    pca(2, pc);

    // 输出前两个主成分
    printf("Principal components:\n");
    for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < D; ++j) {
            printf("%.4lf ", pc[i*D+j]);
        }
        printf("\n");
    }

    return 0;
}

其中，我使用了Eigen库来求解特征值和特征向量，在Linux环境下可以通过以下命令安装：

sudo apt install libeigen3-dev