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

SugarCane88

请教各位高手：

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

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

1. #include <iostream>
   
2. #include "Eigen/Dense"
   
3. #include "pca.h"
   
4. #include <fstream>
   
5. 
   
6. using Eigen::MatrixXd;
   
7. using namespace std;
   
8. 
   
9. int main()
   
10. {
    
11. 
    
12.         int rows = 150; //number of samples
    
13.         int columns = 4; //number of components
    
14.         string pathData = "test_data.txt";
    
15. 
    
16.         MatrixXd m(rows, columns);
    
17. 
    
18.         ifstream istrm;
    
19.         istrm.open(pathData);
    
20. 
    
21.         double a{ 0 };
    
22.         for (int i = 0; i < 150; i++)
    
23.         {
    
24.                 istrm >> a;
    
25.                 m(i, 0) = a;
    
26.                 istrm >> a;
    
27.                 m(i, 1) = a;
    
28.                 istrm >> a;
    
29.                 m(i, 2) = a;
    
30.                 istrm >> a;
    
31.                 m(i, 3) = a;
    
32.         }
    
33.         istrm.close();
    
34. 
    
35. 
    
36. 
    
37.         //for (int i = 0; i < 150; i++)
    
38.         //{
    
39.         //        printf("m的值：%lf,%lf,%lf,%lf\n", m(i, 0), m(i, 1), m(i, 2), m(i, 3));
    
40.         //}
    
41. 
    
42.         MatrixXd uReduced = PCA::Compute(m);
    
43.         MatrixXd mReduced = m*uReduced;
    
44. 
    
45.         cout << uReduced << endl;
    
46.         system("pause");
    
47. }

复制代码

上面是主函数

1. #ifndef PCA_H
   
2. #define PCA_H
   
3. 
   
4. #include <iostream>
   
5. #include <assert.h>
   
6. #include <vector>
   
7. #include "Eigen/Dense"
   
8. #include "Eigen/Eigenvalues"
   
9. 
   
10. using namespace Eigen;
    
11. using namespace std;
    
12. 
    
13. class PCA
    
14. {
    
15. public:
    
16. 
    
17.         /*
    
18.          * Computes the princcipal component of a given matrix. Computation steps:
    
19.          * Compute the mean image.
    
20.          * Compute standard dispersion
    
21.          * Normalizing input data
    
22.          * Compute the covariance matrix
    
23.          * Calculate the eigenvalues and eigenvectors for the covariance matrix
    
24. 
    
25.          * @input MatrixXd D the data samples matrix.
    
26.          *
    
27.          * @returns MatrixXd The eigenVectors of principal component
    
28.          */
    
29. 
    
30.         static MatrixXd Compute(MatrixXd D)
    
31.         {
    
32.                 // 1. Compute the mean image
    
33.                 MatrixXd mean(1, D.cols());
    
34.                 mean.setZero();
    
35. 
    
36.                 for (int i = 0; i < D.cols(); i++)
    
37.                 {
    
38.                         mean(0, i) = D.col(i).mean();
    
39.                 }
    
40. 
    
41. 
    
42.                 // 2. Compute standard dispersion
    
43.                 MatrixXd s(1, D.cols());
    
44.                 for (int i = 0; i < D.cols(); i++)
    
45.                 {
    
46.                         double ss{ 0 };
    
47. 
    
48.                         for (int j = 0; j < D.rows(); j++)
    
49.                         {
    
50.                                 double tmp = (D(j, i) - mean(0, i));
    
51.                                 tmp *= tmp;
    
52.                                 ss += tmp;
    
53.                         }
    
54.                         ss /= (D.rows() - 1);
    
55.                         s(0, i) = sqrt(ss);
    
56.                 }
    
57.                 //cout << s << endl;
    
58. 
    
59. 
    
60.                 // 3. Normalizing input data
    
61.                 MatrixXd normD(D);
    
62.                 normD.setZero();
    
63.                 for (int i = 0; i < D.cols(); i++)
    
64.                 {
    
65.                         for (int j = 0; j < D.rows(); j++)
    
66.                         {
    
67.                                 normD(j, i) = (D(j, i) - mean(0, i)) / s(0, i);
    
68.                         }
    
69.                 }
    
70.                 //cout << normD << endl;
    
71. 
    
72. 
    
73.                 // 3. Compute the covariance matrix
    
74.                 mean.setZero();
    
75. 
    
76.                 for (int i = 0; i < normD.cols(); i++)
    
77.                 {
    
78.                         mean(0, i) = normD.col(i).mean();
    
79.                 }
    
80. 
    
81.                 MatrixXd covariance(D.cols(), D.cols());
    
82.                 covariance.setZero();
    
83.                 for (int k = 0; k < D.cols(); k++)
    
84.                 {
    
85.                         for (int j = 0; j < D.cols(); j++)
    
86.                         {
    
87.                                 double sum{ 0 };
    
88.                                 for (int i = 0; i < D.rows(); i++)
    
89.                                 {
    
90.                                         sum += (normD(i, k) - mean(0, k))*(normD(i, j) - mean(0, j));
    
91.                                 }
    
92.                                 covariance(k, j) = sum / (D.rows() - 1);
    
93.                         }
    
94.                 }
    
95. 
    
96. 
    
97.                 covariance(0, 0) = 0.685693512304251;
    
98.                 covariance(0, 1) = -0.0424340044742729;
    
99.                 covariance(0, 2) = 1.27431543624161;
    
100.                 covariance(0, 3) = 0.516270693512304;
     
101. 
     
102.                 covariance(1, 0) = -0.0424340044742729;
     
103.                 covariance(1, 1) = 0.189979418344519;
     
104.                 covariance(1, 2) = -0.329656375838926;
     
105.                 covariance(1, 3) = -0.121639373601790;
     
106. 
     
107.                 covariance(2, 0) = 1.27431543624161;
     
108.                 covariance(2, 1) = -0.329656375838926;
     
109.                 covariance(2, 2) = 3.11627785234899;
     
110.                 covariance(2, 3) = 1.29560939597315;
     
111. 
     
112.                 covariance(3, 0) = 0.516270693512304;
     
113.                 covariance(3, 1) = -0.121639373601790;
     
114.                 covariance(3, 2) = 1.29560939597315;
     
115.                 covariance(3, 3) = 0.581006263982103;
     
116. 
     
117.                 //cout << "斜方差" << "\n"<<endl;
     
118.                 //cout << covariance << endl;
     
119. 
     
120. 
     
121.                 // 4. Calculate the eigenvalues and eigenvectors for the covariance matrix
     
122.                 EigenSolver<MatrixXd> solver(covariance);
     
123.                 MatrixXd eigenVectors = solver.eigenvectors().real();
     
124.                 VectorXd eigenValues = solver.eigenvalues().real();
     
125. 
     
126.                 cout << eigenVectors << endl;
     
127.                 cout << "\n";
     
128. 
     
129.                 cout << eigenValues << endl;
     
130.                 cout << "\n";
     
131. 
     
132. 
     
133.                 // 5. Find out an number of most important principal components with broken stick method   
     
134.                 int newDimension = brockenStickMethod(eigenValues);
     
135. 
     
136.                 // 6. Return the matrix of eigenVectors appropriate most important components
     
137.                 MatrixXd featureVectors(D.cols(), newDimension);
     
138. 
     
139.                 VectorXd sortEV(eigenValues);
     
140.                 sort(sortEV.derived().data(), sortEV.derived().data() + sortEV.derived().size());
     
141. 
     
142.                 for (int i = 0; i < newDimension; i++)
     
143.                 {
     
144.                         double max = sortEV(sortEV.rows() - 1 - i);
     
145. 
     
146.                         for (int j = 0; j < eigenValues.size(); j++)
     
147.                         {
     
148.                                 if (max == eigenValues(j))
     
149.                                 {
     
150.                                         VectorXd tmp = eigenVectors.col(j);
     
151.                                         featureVectors.col(i) = tmp;
     
152.                                 }
     
153.                         }
     
154.                 }
     
155.                
     
156.                 return featureVectors;
     
157.         }
     
158. 
     
159. private:
     
160.         static int brockenStickMethod(VectorXd eigenValues)
     
161.         {
     
162.                 // Compute number of components, which most influed
     
163.                 double trc{ eigenValues.sum() };
     
164.                 VectorXd l(eigenValues.size());
     
165.                 int countPrincComp{ 0 };
     
166.                 l.setZero();
     
167. 
     
168.                 for (int i = 0; i < eigenValues.size(); i++)
     
169.                 {
     
170.                         double sum{ 0 };
     
171.                         for (int j = i; j < eigenValues.size(); j++)
     
172.                         {
     
173.                                 double coef = 1;
     
174.                                 sum += (coef / j);
     
175.                         }
     
176.                         l(i) = sum / eigenValues.size();
     
177.                 }
     
178. 
     
179. 
     
180.                 VectorXd sortEV(eigenValues);
     
181.                 sort(sortEV.derived().data(), sortEV.derived().data() + sortEV.derived().size());
     
182. 
     
183.                 for (int i = 0; i < sortEV.size(); i++)
     
184.                 {
     
185.                         if ((sortEV(i) / trc) >l(i)) { countPrincComp++; }
     
186.                 }
     
187. 
     
188.                 return countPrincComp;
     
189.         }
     
190. };
     
191. 
     
192. #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等软件算出的特征值一致。