C++版本主成分分析（PCA）的特征向量与Matlab自带PCA函数特征向量正负号不同，求助。

SugarCane88

请教各位高手：

        用C++编写的主成分分析（PCA）求解得到的特征向量与Matlab自带PCA函数特征向量正负号有些不同。
        特征向量数值是相同的，特征值也一样，就特征向量的正负符号有些是相反的。
        我在网上找过两个自编版本，都出现类似问题。

        困扰好久了，请高人指点。

#include <iostream>
#include "Eigen/Dense"
#include "pca.h"
#include <fstream>

using Eigen::MatrixXd;
using namespace std;

int main()
{

        int rows = 150; //number of samples 
        int columns = 4; //number of components
        string pathData = "test_data.txt";

        MatrixXd m(rows, columns);

        ifstream istrm;
        istrm.open(pathData);

        double a{ 0 };
        for (int i = 0; i < 150; i++)
        {
                istrm >> a;
                m(i, 0) = a;
                istrm >> a;
                m(i, 1) = a;
                istrm >> a;
                m(i, 2) = a;
                istrm >> a;
                m(i, 3) = a;
        }
        istrm.close();



        //for (int i = 0; i < 150; i++)
        //{
        //        printf("m的值：%lf,%lf,%lf,%lf\n", m(i, 0), m(i, 1), m(i, 2), m(i, 3));
        //}

        MatrixXd uReduced = PCA::Compute(m);
        MatrixXd mReduced = m*uReduced;

        cout << uReduced << endl;
        system("pause");
}

上面是主函数

#ifndef PCA_H
#define PCA_H

#include <iostream>
#include <assert.h>
#include <vector>
#include "Eigen/Dense"
#include "Eigen/Eigenvalues"

using namespace Eigen;
using namespace std;

class PCA
{
public:

        /*
         * Computes the princcipal component of a given matrix. Computation steps:
         * Compute the mean image.
         * Compute standard dispersion
         * Normalizing input data
         * Compute the covariance matrix
         * Calculate the eigenvalues and eigenvectors for the covariance matrix

         * @input MatrixXd D the data samples matrix.
         *
         * @returns MatrixXd The eigenVectors of principal component 
         */

        static MatrixXd Compute(MatrixXd D)
        {
                // 1. Compute the mean image 
                MatrixXd mean(1, D.cols());
                mean.setZero();

                for (int i = 0; i < D.cols(); i++)
                {
                        mean(0, i) = D.col(i).mean();
                }


                // 2. Compute standard dispersion
                MatrixXd s(1, D.cols());
                for (int i = 0; i < D.cols(); i++)
                {
                        double ss{ 0 };

                        for (int j = 0; j < D.rows(); j++)
                        {
                                double tmp = (D(j, i) - mean(0, i));
                                tmp *= tmp;
                                ss += tmp;
                        }
                        ss /= (D.rows() - 1);
                        s(0, i) = sqrt(ss);
                }
                //cout << s << endl;


                // 3. Normalizing input data 
                MatrixXd normD(D);
                normD.setZero();
                for (int i = 0; i < D.cols(); i++)
                {
                        for (int j = 0; j < D.rows(); j++)
                        {
                                normD(j, i) = (D(j, i) - mean(0, i)) / s(0, i);
                        }
                }
                //cout << normD << endl;


                // 3. Compute the covariance matrix
                mean.setZero();

                for (int i = 0; i < normD.cols(); i++)
                {
                        mean(0, i) = normD.col(i).mean();
                }

                MatrixXd covariance(D.cols(), D.cols());
                covariance.setZero();
                for (int k = 0; k < D.cols(); k++)
                {
                        for (int j = 0; j < D.cols(); j++)
                        {
                                double sum{ 0 };
                                for (int i = 0; i < D.rows(); i++)
                                {
                                        sum += (normD(i, k) - mean(0, k))*(normD(i, j) - mean(0, j));
                                }
                                covariance(k, j) = sum / (D.rows() - 1);
                        }
                }


                covariance(0, 0) = 0.685693512304251;
                covariance(0, 1) = -0.0424340044742729;
                covariance(0, 2) = 1.27431543624161;
                covariance(0, 3) = 0.516270693512304;

                covariance(1, 0) = -0.0424340044742729;
                covariance(1, 1) = 0.189979418344519;
                covariance(1, 2) = -0.329656375838926;
                covariance(1, 3) = -0.121639373601790;

                covariance(2, 0) = 1.27431543624161;
                covariance(2, 1) = -0.329656375838926;
                covariance(2, 2) = 3.11627785234899;
                covariance(2, 3) = 1.29560939597315;

                covariance(3, 0) = 0.516270693512304;
                covariance(3, 1) = -0.121639373601790;
                covariance(3, 2) = 1.29560939597315;
                covariance(3, 3) = 0.581006263982103;

                //cout << "斜方差" << "\n"<<endl;
                //cout << covariance << endl;


                // 4. Calculate the eigenvalues and eigenvectors for the covariance matrix
                EigenSolver<MatrixXd> solver(covariance);
                MatrixXd eigenVectors = solver.eigenvectors().real();
                VectorXd eigenValues = solver.eigenvalues().real();

                cout << eigenVectors << endl;
                cout << "\n";

                cout << eigenValues << endl;
                cout << "\n";


                // 5. Find out an number of most important principal components with broken stick method   
                int newDimension = brockenStickMethod(eigenValues);

                // 6. Return the matrix of eigenVectors appropriate most important components
                MatrixXd featureVectors(D.cols(), newDimension);

                VectorXd sortEV(eigenValues);
                sort(sortEV.derived().data(), sortEV.derived().data() + sortEV.derived().size());

                for (int i = 0; i < newDimension; i++)
                {
                        double max = sortEV(sortEV.rows() - 1 - i);

                        for (int j = 0; j < eigenValues.size(); j++)
                        {
                                if (max == eigenValues(j))
                                {
                                        VectorXd tmp = eigenVectors.col(j);
                                        featureVectors.col(i) = tmp;
                                }
                        }
                }
                
                return featureVectors;
        }

private:
        static int brockenStickMethod(VectorXd eigenValues)
        {
                // Compute number of components, which most influed
                double trc{ eigenValues.sum() };
                VectorXd l(eigenValues.size());
                int countPrincComp{ 0 };
                l.setZero();

                for (int i = 0; i < eigenValues.size(); i++)
                {
                        double sum{ 0 };
                        for (int j = i; j < eigenValues.size(); j++)
                        {
                                double coef = 1;
                                sum += (coef / j);
                        }
                        l(i) = sum / eigenValues.size();
                }


                VectorXd sortEV(eigenValues);
                sort(sortEV.derived().data(), sortEV.derived().data() + sortEV.derived().size());

                for (int i = 0; i < sortEV.size(); i++)
                {
                        if ((sortEV(i) / trc) >l(i)) { countPrincComp++; }
                }

                return countPrincComp;
        }
};

#endif // PCA_H

这个是pac.h文件。

在测试阶段，我把主成分步骤中corr的值替换成cov，使得它和Matlab版本一致。

SugarCane88

我回答一下我自己的帖子，实在是花了太多的时间：

并不是编程的问题，而是：

主成分的特征向量百有两个约束条件：（1）特征向量的模为1；（2）特征向量度两两正交。在这两个条件的制约下，一问个特征值对应两个方向相反的特征向量a和-a。因此答需要再设定一个约束条件，即：取值最大的样本的主成回分的得分必须大于取值最小的样本的主成分的得分，满足这个条件的特答征向量就只有一个了。

SugarCane88

SugarCane88 发表于 2020-5-6 00:37
我回答一下我自己的帖子，实在是花了太多的时间：

并不是编程的问题，而是：

为了给后续寻找这个问题答案的人少走弯路，我再次回答：

主要就是这句话“forcing the largest component of each vector to be positive”
即如果想要 sign of egienvector唯一的话，确保各个成分中最大的值为正，如某个主成分向量中绝对值最大的值为负，则整行取相反即可。

例如：
主成分系数为“0.603033336168387
-0.735701512298084
-0.281699372506027
-0.125457338572313”

由于绝对值最大的“-0.735701512298084”为负，所以该主成分对应的特征向量全部取反：
“-0.603033336168387
0.735701512298084
0.281699372506027
0.125457338572313”

这样可保证与matlab,python等软件算出的特征值一致。