Python与C++中SVD（奇异值分解）得到的右奇异值不同

目录

一、现象与问题

​二、解决办法

1.C++与Python中SVD的结果不同问题。

2.特征值求解右奇异值（eig）与直接求解右奇异值（svd）结果不同

i.发现有几列元素互为相反数

ii.最后两列的位置互换

iii.最好直接使用奇异值分解，而不是使用分步计算。其中会有很多误差。

---

一、现象与问题

在Python中，使用的是np.linalg.svd()来得到左右奇异值和特征值。

[U,S,V] = np.linalg.svd(A)

其中A是8*9的矩阵，得到的U是左奇异值是8*8的，S是8个值，V是9*9的。

Python中A的矩阵值如下：

其中右奇异值V的数据显示如下：

在C++中，使用的Eigen库中方法

cout<<"AAAAAAAAAAAA"<<endl;
cout<<A<<endl;
Eigen::JacobiSVD<Eigen::MatrixXf> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
cout<<"+++++++++++++++++++++++++"<<endl;
cout<<svd.matrixV()<<endl;

得到右奇异值V却是9*8列的，很是奇怪。

为了得到正确的右奇异值

使用原矩阵ATA进行计算特征向量得到右奇异值，详细原理见奇异值分解

PtshV ATA,D,V,Q;
ATA=A.transpose()*A;
Eigen::EigenSolver<PtshV> eigenSolver(ATA);
D=eigenSolver.pseudoEigenvalueMatrix();
Q=eigenSolver.pseudoEigenvectors();
V=Q.transpose();
cout << "The  D is:" << endl << D << endl;
cout << "The  matrix Q is:" << endl << Q << endl;
cout << "The  matrix Q.transpose is:" << endl << Q.transpose() << endl;

得到的结果：

其中得到的Q.transpose()是Q的转置，即所求的右奇异值，是9*9列的。

但是仔细对比，还是发现使用特征值分解得到的右奇异值与Python中得到的右奇异值不同。

再在Python中使用特征向量计算的方法求得右奇异值

m=np.dot(A.T,A)
a1,a2=np.linalg.eig(m)

这里a2是求得的特征向量的矩阵。

二、解决办法

1.C++与Python中SVD的结果不同问题。

后来在仔细对比C++Eigen中SVD和np.linalg.svd()的源代码，发现了其中的差异。

np.linalg.svd()：

def svd(a, full_matrices=True, compute_uv=True):

其中出来要进行svd处理的矩阵，还有两个bool参数是默认为True的。

full_matrices是用来判断是否要进行全矩阵计算。

full_matrices=True  表明 若a为m*n的矩阵，则u为m*m的矩阵，v为n*n的矩阵。

反之，u为m*k的矩阵，v为k×n的矩阵；其中k为min（m，n）。

compute_uv是用来判断出来S以外是否计算u和v的值。Ture为计算，False不计算。

之前的计算都是取默认值

[U,S,V] = np.linalg.svd(A)

还有一点需要注意：这里的V是特征向量矩阵转置过的，即一行是一个特征向量。

JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)

第一个参数为要进行计算的矩阵，第二个参数有四个取值：#ComputeFullU, #ComputeThinU,#ComputeFullV, #ComputeThinV.

Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf).

只有当矩阵类型具有动态的列时，才能使用精简单位。

将之前的计算代码修改为：

typedef Eigen::Matrix<float ,9,9> PtshV;
PtshV V；

Eigen::JacobiSVD<Eigen::MatrixXf> svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV);
V=svd.matrixV();
cout<<"VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV"<<endl;
cout<<V.transpose()<<endl;

这里的V是特征向量矩阵未进行转置的，一列是一个特征向量。

结果如下：

与之前的Python版本对比，误差忽略不计。

2.特征值求解右奇异值（eig）与直接求解右奇异值（svd）结果不同

eig求得的是特征向量的矩阵，而svd得到的是特征向量矩阵的转置。

将右奇异值V转置a3与eig得到的特征向量的矩阵a2进行比较。

i.发现有几列元素互为相反数

这种现象是正常的，因为这些列向量的本质就是特征向量，同一特征值的特征向量的非零线性组合任然是特征向量，所以-1倍没有错。

ii.最后两列的位置互换

这是因为其特征值的问题，对于奇异值分解的特征值是对两个解法求得的特征值开平方的得到的。详细见此。

但是直接求解特征值得到的右奇异值的方法，得到的第8个特征向量是负数，a1是直接求解到特征值

下图是奇异值分解得到的奇异值的平方，也是特征值。

所以在奇异值分解的时候，若不是full_matrices的求解，就会舍弃第七个负值，而full_matrices=True的求解方法会将这个负值的

特征向量放在最后。

iii.最好直接使用奇异值分解，而不是使用分步计算。其中会有很多误差。

iv.在a2中6行7列和7行7列的值也不同，还未解决。

---

可参考：

https://zhuanlan.zhihu.com/p/43578482

https://codeday.me/bug/20181026/336042.html