PCA降维的C语言实现

已于 2024-03-25 15:58:25 修改
阅读量163 收藏
点赞数 2
文章标签： c语言
于 2024-03-25 15:57:23 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_53383166/article/details/137016966
版权 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int max_row_order(int n, double A[n][n])
{
    double m = 0;
    int mro = 0;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (j == i)
            {
                continue;
            }
            
            if (fabs(A[i][j]) > m)
            {
                m = fabs(A[i][j]);
                mro = i;
            }  
        }                 
    }
    int p = mro;
    return p;
}

int max_column_order(int n, double A[n][n])
{
    double m = 0;
    int mco = 0;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (j == i)
            {
                continue;
            }
            
            if (fabs(A[i][j]) > m)
            {
                m = fabs(A[i][j]);
                mco = j;
            }  
        }                 
    }
    int q = mco;
    return q;
}

double t_slove(int n, double A[n][n], int p, int q)
{
    double s = (A[q][q] - A[p][p])/(2*A[p][q]);
    double t;
    if (s == 0)
    {
        t = 1;
    }
    else
    {
        double x1 = -s + sqrt(s*s+1);
        double x2 = -s - sqrt(s*s+1);
        if (fabs(x1) > fabs(x2))
        {
            t = x2;
        }
        else
        {
            t = x1;
        }
    }
    return t;
}

double calculate_nondiagonal(int n, double A[n][n])
{
    double sum = 0;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (j == i)
            {
                continue;
            }
            sum = sum + A[i][j] * A[i][j];
        }                 
    }
    return sum;
}

void print_matrix(int row, int col, double matrix[row][col])
{
    for(int i = 0;i < row;i ++)
    {
        for(int j = 0;j < col;j ++)
        {
            printf("%.12lf\t",matrix[i][j]);
        }    
        printf("\n");
    }
    printf("----------------------------------------------");
    printf("\n");
}

void print_vector(int n, double vector[n])
{
    for(int i = 0;i < n;i ++)
    {
        printf("%.12lf\n",vector[i]);   
    }
    printf("----------------------------------------------");
    printf("\n");
}

double calculate_mode(int row, double x[row][1])
{
    double mode2 = 0;
    for (int i = 0; i < row; i++)
    {
        mode2 = mode2 + x[i][0]*x[i][0];
    }
    double mode = sqrt(mode2);
    return mode;
}

double * GaussEliminationWithPartialPivoting1(int n, double A[n-1][n-1], double b[n-1])
{
    double *result;
    int size = n-1;
    result = (double *)malloc(size*sizeof(double));     
    for (int i = 0; i < n-1; i++)
    {
        result[i] = 0;
    }

    for (int i = 0; i < n-1; i++)
    {
        int k = i;
        for (int j = i+1; j < n-1; j++)
        {
            if (fabs(A[k][i]) < fabs(A[j][i]))
            {
                k = j;
            }
        }
        for (int m = 0; m < n-1; m++)
        {
            double c;
            c = A[i][m];
            A[i][m] = A[k][m];
            A[k][m] = c;
        }
        double d;
        d = b[i];
        b[i] = b[k];
        b[k] = d;

        for (int j = i+1; j < n-1; j++)
        {
            double v = -A[j][i]/A[i][i];
            for (int k = 0; k < n-1; k++)
            {
                A[j][k] = v*A[i][k] + A[j][k];
            }
            b[j] = v*b[i] + b[j];
        }  
    }
    
    for (int i = n-2; i >= 0; i--)
    {
        if (i == n-2)
        {
            result[i] = 1.0;
        }
        else
        {
            for (int j = n-2; j >= i+1; j--)
            {
                b[i] = b[i] - A[i][j]*result[j];
            }
            result[i] = b[i]/A[i][i];
        }    
    }
    printf("After Elimination \n");
    for(int i = 0;i < n-1;i ++)
    {
        for(int j = 0;j < n-1;j ++)
        {
            printf("%lf\t",A[i][j]);
        }    
        printf("\n");
    }
    printf("----------------------------------------------");
    printf("\n");
    return result;
    free(result);
}

double calculate_mean (int m, double vector[m])
{
    double sum = 0;
    for (int i = 0; i < m; i++)
    {
        sum = sum + vector[i];
    }
    double mean = sum/m;
    return mean;
}

int main()
{
    FILE *fp = NULL;
    fp = fopen("iris.txt", "r");
    int m = 150;
    int n = 4;
    double a[m], b[m], c[m], d[m], e[m];
    for (int i = 0; i < m; i++)
    {
        fscanf(fp, "%lf,%lf,%lf,%lf,%lf", &a[i], &b[i],&c[i],&d[i],&e[i]);
    }
    double mean_a = calculate_mean(m, a);
    double mean_b = calculate_mean(m, b);
    double mean_c = calculate_mean(m, c);
    double mean_d = calculate_mean(m, d);
    double X[n][m];
    for (int i = 0; i < n; i++)//decentralize X
    {
        for (int j = 0; j < m; j++)
        {
            if (i == 0)
            {
                X[i][j] = a[j] - mean_a;
            }
            else if (i == 1)
            {
                X[i][j] = b[j] - mean_b;
            }
            else if (i == 2)
            {
                X[i][j] = c[j] - mean_c;
            }
            else if (i == 3)
            {
                X[i][j] = d[j] - mean_d;
            }
            
        }
        
    }
  
    //calculate XXT/m -> A
    double A[n][n];
    for (int i = 0; i < n; i++)//A = 0
    {
        for (int j = 0; j < n; j++)
        {
            A[i][j] = 0;
        }   
    }
    for (int i = 0; i < n; i++)//A = XXT/m
    {
        for (int j = 0; j < n; j++)
        {
            for (int k = 0; k < m; k++)
            {
                A[i][j] = A[i][j] + X[i][k]*X[j][k]/m;
            }            
        }   
    }
    double Ao[n][n];
    for (int i = 0; i < n; i++)//Ao = A
    {
        for (int j = 0; j < n; j++)
        {            
            Ao[i][j] = A[i][j];                      
        }   
    }
    printf("A (XXT/m) = \n");
    print_matrix(n, n, A);

    printf("=======================================================================\n");
    //calculate eigenvalue of A
    double limit = pow(10, -6);
    int p, q;
    double t, cos, sin;
    double Apio, Aqio;

    printf("calculate eigenvalue of A\n");
    int num = 0;
    while (calculate_nondiagonal(n, A) > limit)
    {        
        printf("sum of square of nondiagonal elements = %.12lf\n", calculate_nondiagonal(n, A));
        p = max_row_order(n,A);
        q = max_column_order(n,A);
        t = t_slove(n,A, p, q);
        cos = 1/sqrt(1 + t*t);
        sin = t/sqrt(1 + t*t);
        printf("Try %d \n", num+1);
        printf("choose p =%d, q = %d \n", p+1, q+1);
        printf("t = %lf, cos = %lf, sin = %lf \n", t, cos, sin);
        
        A[p][p] = A[p][p] - t*A[p][q];
        A[q][q] = A[q][q] + t*A[p][q];
        A[p][q] = 0;
        A[q][p] = 0;
        for (int i = 0; i < n; i++)
        {            
            if (i != p && i!= q)
            {
                Apio = A[p][i];
                Aqio = A[q][i];
                A[i][p] = cos*Apio - sin*Aqio;
                A[p][i] = cos*Apio - sin*Aqio;
                A[i][q] = sin*Apio + cos*Aqio;
                A[q][i] = sin*Apio + cos*Aqio;
            }                       
        }
        print_matrix(n, n, A);
        num ++;
    }
    double eigvalue[n][1];
    for (int i = 0; i < n; i++)
    {
        eigvalue[i][0] = A[i][i];
    }
    //sort eigenvalue
    for(int i = 0; i < n-1; i++)
    {
        int m = i;
        for(int j = i+1; j < n; j++)
        {
            if(eigvalue[j][0] > eigvalue[m][0])
                m = j;
        }            
        if(m != i)
        {
            double temp = eigvalue[i][0];
            eigvalue[i][0] = eigvalue[m][0];
            eigvalue[m][0] = temp;
        }
    }
    printf("eigenvalue of A = \n");
    print_matrix(n, 1, eigvalue);
    
    printf("=======================================================================\n");
    printf("choose the two biggest eigenvalue, and calculate the corresponding unit eigenvectors:\n");   
    //calculte unit eigenvector of A
    double I[n][n];
    for (int i = 0; i < n; i++)//initialize I
    {
        for (int j = 0; j < n; j++)
        {            
            if (i == j)
            {
                I[i][j] = 1;
            }
            else
            {
                I[i][j] = 0;
            }             
        }   
    }
    double bo[n];   
    double U[n][2];
    for (int ia = 0; ia < 2; ia++)
    {
        printf("when λ = %.12lf\n", eigvalue[ia][0]);
        double B[n][n];
        for (int i = 0; i < n; i++)//B = A - λΙ
        {
            for (int j = 0; j < n; j++)
            {
                B[i][j] = Ao[i][j] - I[i][j]*eigvalue[ia][0];
            }   
        }
        printf("A - λΙ = \n");
        print_matrix(n, n, B);    

        for (int i = 0; i < n; i++)//b = 0
        {
            bo[i] = 0;
        }
        double *result = GaussEliminationWithPartialPivoting1(n+1, B, bo);
        double u[n][1];
        for (int i = 0; i < n; i++)//u = result
        {
            u[i][0] = result[i];
        } 
        double mode = calculate_mode(n, u);       
        for (int i = 0; i < n; i++)//u = result
        {
            u[i][0] = u[i][0]/mode;
        } 
        printf("corresponding eigenvector = \n");
        print_matrix(n, 1, u);  
        for (int i = 0; i < n; i++)//U
        {
            U[i][ia] = u[i][0];               
        }
    }
    printf("U, composed by eigenvectors of AAT, = \n");
    print_matrix(n, 2, U);     
    
    /*==================================================================================*/
    printf("=======================================================================\n");
    //calculate projection
    double proj1[m];
    for (int i = 0; i < m; i++)//projection = 0
    {
        proj1[i] = 0;        
    }
    for (int i = 0; i < m; i++)//projection on eigenvector1
    {
        for (int j = 0; j < n; j++)
        {
            proj1[i] = proj1[i] + U[j][0]*X[j][i];
        }        
    }
    double proj2[m];
    for (int i = 0; i < m; i++)//projection = 0
    {
        proj2[i] = 0;        
    }
    for (int i = 0; i < m; i++)//projection on eigenvector2
    {
        for (int j = 0; j < n; j++)
        {
            proj2[i] = proj2[i] + U[j][1]*X[j][i];
        }        
    }
    printf("projection1: \n");
    print_vector(m, proj1);
    printf("projection2: \n");
    print_vector(m, proj2);
    //
    printf("these data would be presented RED\n");
    printf("projection1\tprojection2\n");
    FILE *fr = NULL;
    fr = fopen("RED.txt", "w+");
    for (int i = 0; i < m; i++)
    {
        if (e[i] == 0)
        {
            printf("%lf\t%lf\n", proj1[i], proj2[i]);
            fprintf(fr, "%lf %lf\n", proj1[i], proj2[i]);
        }        
    }
    fclose(fr);
    //
    printf("these data would be presented GREEN\n");
    printf("projection1\tprojection2\n");
    FILE *fg = NULL;
    fg = fopen("GREEN.txt", "w+");
    for (int i = 1; i < m; i++)
    {
        if (e[i] == 1)
        {
            printf("%lf\t%lf\n", proj1[i], proj2[i]);
            fprintf(fg, "%lf %lf\n", proj1[i], proj2[i]);
        }        
    }
    fclose(fg);
    //
    printf("these data would be presented BLUE\n");
    printf("projection1\tprojection2\n");
    FILE *fb = NULL;
    fb = fopen("BLUE.txt", "w+");
    for (int i = 0; i < m; i++)
    {
        if (e[i] == 2)
        {
            printf("%lf\t%lf\n", proj1[i], proj2[i]);
            fprintf(fb, "%lf %lf\n", proj1[i], proj2[i]);
        }        
    }
    fclose(fb);
      
    fclose(fp);
    return 0;
}
作业不能白做
关注
  • 2
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
pca算法 c语言,降维算法一 : PCA (Principal Component Analysis)
weixin_33293858的博客
05-25 763
PCA简介在许多领域的研究与应用中，往往需要对反映事物的多个变量进行大量的观测，收集大量数据以便进行分析寻找规律。多变量大样本无疑会为研究和应用提供了丰富的信息，但也在一定程度上增加了数据采集的工作量，更重要的是在多数情况下，许多变量之间可能存在相关性，从而增加了问题分析的复杂性，同时对分析带来不便。如果分别对每个指标进行分析，分析往往是孤立的，而不是综合的。盲目减少指标会损失很多信息，容易产生错...
C/C++实现PCA降维以及故障监测
Oldwater1的博客
10-28 1508
C/C++实现PCA降维以及故障监测@TOC 欢迎使用Markdown编辑器 你好！ 这是你第一次使用 Markdown编辑器 所展示的欢迎页。如果你想学习如何使用Markdown编辑器, 可以仔细阅读这篇文章，了解一下Markdown的基本语法知识。 新的改变 我们对Markdown编辑器进行了一些功能拓展与语法支持，除了标准的Markdown编辑器功能，我们增加了如下几点新功能，帮助你用它写博客： 全新的界面设计 ，将会带来全新的写作体验； 在创作中心设置你喜爱的代码高亮样式，Markdown 将代码
PCA算法C语言描述
07-16
C语言实现的主成分分析法，PCA算法用于各个领域妥妥的
PCA的C语言代码
boyhailong的专栏
12-02 7414
来自维基百科： C*********************** Contents **************************************** C* Principal Components Analysis: C, 638 lines. **************************** C* Sample input data set (final 36 li
PCA降维及C代码
EatAppleS的专栏
03-03 1633
转载至：http://blog.csdn.net/u011001084/article/details/51363892转载至：http://blog.csdn.net/jinshengtao/article/details/18599165PCA人脸识别        将PCA用于人脸识别的过程如下：        1.假设有400幅尺寸为100*100的图像，构成10000*400的矩阵；  ...
pca—SIFT C语言
11-18
PCA—SIFT demo 主要是有C的文件为来match。可以使用的、
pca算法 c语言,Opencv中的pca算法
weixin_35664144的博客
05-25 408
对于PCA，一直都是有个概念，没有实际使用过，今天终于实际使用了一把，发现PCA还是挺神奇的。在OPENCV中使用PCA非常简单，只要几条语句就可以了。1、初始化数据//每一行表示一个样本CvMat* pData = cvCreateMat( 总的样本数, 每个样本的维数, CV_32FC1 );CvMat* pMean = cvCreateMat(1, 样本的维数, CV_32FC1);//pE...
PCA降维处理C++代码
10-20
此程序是一个简单的PCA降维处理，程序中以SampleData.txt中的三维数据点作为输入，程序的输出存放在Results目录， 其中包括： 1. PCAResult.txt 存放经主元分析后在主元轴坐标系下的数据坐标。 2. DisplaySamples.ms 为3ds Max脚本文件，用3ds Max打开后可看到SampleData.txt表示的数据点。 3. DisplayResult.ms 为3ds Max脚本文件，用3ds Max打开后可看到将原3维数据在主元轴坐标系下降到2维后的数据点；同时该文件还显示了主元分析得到的前两个主元轴。
PCA算法cpp实现
windroid.xyz
05-24 414
PCA过程： 零均值化 求协方差 求特征值及特征向量 按特征值大小排序特征向量 取前k行组成变换矩阵 使用变换矩阵即可进行降维或还原 测试数据集（68040*32）：CorelFeatures-mld/ColorHistogram.asc #include &lt;iostream&gt; #include &lt;math.h&gt; #include &lt;string.h&gt;...
PCA.rar_PCA 特征降维_PCA降维 C_PCA降维C_特征向量 c_特征降维
09-19
在给定的"PCA.rar"压缩包中，可能包含了一些关于如何使用C语言实现PCA降维的代码或者教程。"PCA2"可能是其中的一个子文件，可能包含了实际的PCA算法实现，或者是经过PCA处理后的低维特征向量数据。使用C语言实现PCA...
PCA代码汇总
05-29
因为修改算法需要，找到了一些PCA代码，上传上来，希望对大家有帮助
pca模块c代码
03-30
pca模块c代码，基于C8051f360单片机
PCA代码实现详解
06-14
对主成分分析（PCA）的C++代码实现，里面有对PCA实现步骤的详细讲解，并把自己写的PCA代码与Opencv自带的PCA函数的运行结果进行了对比。
PCA算法C++实现
07-31
使用雅克比算法实现PCA降维
PCA_sourcecode.zip_C 降维_PCA 降维_pca降维_特征降维_降维C
09-23
在C语言实现PCA的过程中，通常包括以下几个关键步骤： 1. **数据预处理**：首先，对原始数据进行标准化或归一化，确保所有特征在同一尺度上。这是因为PCA的结果受数据量纲影响，统一尺度可以避免量纲不同的特征对...
eigen PCA降维算法
01-17
PCA（主成分分析，Principal Component Analysis）是一种广泛应用的数据降维技术，它通过线性变换将原始数据转换到一个新的坐标系中，使得新的坐标系中的各轴是按照数据方差从大到小排列的，从而保留了数据的主要...
PCA人脸识别学习及C语言实现
热门推荐
Where there is life, there is hope
01-21 3万+
人脸识别主要方法：   .Eigenfaces，PCA（Principal Component Analysis），Turk and Pentland，1991   .Fisherfaces，LDA（Linear Discriminant Analysis），Belhumeur, Hespanha and Kriegman，1997   .LBPH，Local Binary Patter
PCA,LDA的C++实现（内附代码~）
ZSZ_shsf的博客
03-31 3038
在一个项目中遇到了这些问题：Visual C++ 无法解析的外部符号问题 然后调用SVD，LDA pca的C++实现
PCA降维matlab实现
最新发布
09-15
在MATLAB中，你可以使用内置的函数`pca`来实现PCA降维。下面是一个简单的示例代码： ```matlab % 假设你的数据存储在一个矩阵X中，每行表示一个样本，每列表示一个特征 X = [1, 2, 3; 4, 5, 6; 7, 8, 9]; % 执行PCA降维操作 [coeff, score, latent] = pca(X); % coeff是主成分（特征向量）矩阵，每列对应一个主成分 % score是降维后的数据矩阵，每行对应一个样本，每列对应一个主成分的得分 % latent是特征值向量，表示每个主成分的方差解释程度 % 假设你想降维到2维 k = 2; reduced_X = X * coeff(:, 1:k); % 输出降维后的数据矩阵 disp(reduced_X); ``` 在这个示例中，`X`是一个3x3的矩阵，我们使用`pca`函数对其进行PCA降维。通过调整`k`的值，你可以选择降维到的维度。 希望这个示例对你有所帮助！如果你有其他问题，请随时提问。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

作业不能白做

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值