numpy opencv matlab eigen SVD结果对比

阅读量：6295 次

发布时间：2019-06-22

本文共 5037 字，大约阅读时间需要 16 分钟。

参考

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

https://byjiang.com/2017/11/18/SVD/

http://www.bluebit.gr/matrix-calculator/

https://stackoverflow.com/questions/3856072/single-value-decomposition-implementation-c

https://stackoverflow.com/questions/35665090/svd-matlab-implementation

矩阵奇异值分解简介及C++/OpenCV/Eigen的三种实现

https://blog.csdn.net/fengbingchun/article/details/72853757

numpy.linalg.svd 源码

https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.html

计算矩阵的特征值与特征向量的方法

https://blog.csdn.net/Junerror/article/details/80222540

https://jingyan.baidu.com/article/27fa7326afb4c146f8271ff3.html

不同的库计算结果不一致

原因在于特征向量不唯一，特征值是唯一的

来源

https://stackoverflow.com/questions/35665090/svd-matlab-implementation

Both are correct... The rows of the v you got from numpy are the eigenvectors of M.dot(M.T) (the transpose would be a conjugate transpose in the complex case). Eigenvectors are in the general case defined only up to a multiplicative constant, so you could multiply any row of v by a different number, and it will still be an eigenvector matrix.

There is the additional constraint on v that it be a , which loosely translates to its rows being orthonormal. This reduces your available choices for every eigenvector to only 2: the normalized eigenvector pointing in either direction. But you still get to multiply any row by -1 and still have a valid v.

A = U * S * V

1 手动计算

给定一个大小为 $2\times 2$ 的矩阵 $A=\left[ \begin{array}{cc} 4 & 4 \\ -3 & 3 \\ \end{array} \right]$ ，虽然这个矩阵是方阵，但却不是对称矩阵，我们来看看它的奇异值分解是怎样的。

由 $AA^T=\left[ \begin{array}{cc} 32 & 0 \\ 0 & 18 \\ \end{array} \right]$ 进行对称对角化分解，得到特征值为 $\lambda_1=32$ ， $\lambda_2=18$ ，相应地，特征向量为 $\vec{p}_1=\left( 1,0 \right) ^T$ ， $\vec{p}_2=\left(0,1\right)^T$ ；由 $A^TA=\left[ \begin{array}{cc} 25 & 7 \\ 7 & 25 \\ \end{array} \right]$ 进行对称对角化分解，得到特征值为 $\lambda_1=32$ ， $\lambda_2=18$ ，

当 lamda1 = 32

AA.T - lamda1*E = [-7 7

　　　　　　　　　7 -7]

线性变换【-1 1

　　　　　0 0】

-x1 + x2 = 0

x1 = 1 x2 = 1 特征向量为【1 1】.T 归一化为【1/sqrt(2), 1/sqrt(2)】

x1 = -1 x2 = -1 特征向量为【-1 -1】.T 归一化为【-1/sqrt(2), -1/sqrt(2)】

当 lamda2 = 18

AA.T - lamda2*E = [7 7

　　　　　　　　　7 7]

线性变换【1 1

　　　　　0 0】

x1 + x2 = 0

如果x1 = -1 x2 = 1 特征向量为【-1 1】.T 归一化为【-1/sqrt(2), 1/sqrt(2)】

如果 x1 = 1 x2 = -1 特征向量为【-1 1】.T 归一化为【1/sqrt(2), -1/sqrt(2)】可见特征向量不唯一

特征向量为 $\vec{q}_1=\left(\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2}\right)^T$ ， $\vec{q}_2=\left(-\frac{\sqrt{2}}{2}, \frac{\sqrt{2}}{2}\right)^T$ 。取 $\Sigma =\left[ \begin{array}{cc} 4\sqrt{2} & 0 \\ 0 & 3\sqrt{2} \\ \end{array} \right]$ ，则矩阵 $A$ 的奇异值分解为

$A=P\Sigma Q^T=\left(\vec{p}_1,\vec{p}_2\right)\Sigma \left(\vec{q}_1,\vec{q}_2\right)^T$

$=\left[ \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array} \right] \left[ \begin{array}{cc} 4\sqrt{2} & 0 \\ 0 & 3\sqrt{2} \\ \end{array} \right] \left[ \begin{array}{cc} \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} \\ -\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} \\ \end{array} \right] =\left[ \begin{array}{cc} 4 & 4 \\ -3 & 3 \\ \end{array} \right]$ .

2 MATLAB 结果与手动计算不同

AB =  [[ 4 4 ],  [-3   3 ]][U,S,V] = svd(AB);USV'#需要转置

AB =

4 4

-3 3

U =

-1 0

0 1

S =

5.6569 0

0 4.2426

V =

-0.7071 -0.7071

-0.7071 0.7071

3 NUMPY结果与手动计算不同, 与matlab相同它们都是调用lapack的svd分解算法。

import numpy as npimport mathimport cv2a = np.array([[4,4],[-3,3]])# a = np.random.rand(2,2) * 10print(a)u, d, v = np.linalg.svd(a)print(u)print(d)print(v)#不需要转置

A = [[ 4 4]

[-3 3]]

U =

[[-1. 0.]

[ 0. 1.]]

[5.65685425 4.24264069]

[[-0.70710678 -0.70710678]

[-0.70710678 0.70710678]]

4 opencv结果与手动计算结果相同

import numpy as npimport mathimport cv2a = np.array([[4,4],[-3,3]],dtype=np.float32)# a = np.random.rand(2,2) * 10print(a)S1, U1, V1 = cv2.SVDecomp(a)print(U1)print(S1)print(V1)#不需要转置

a = [[ 4. 4.]

[-3. 3.]]

U =

[[0.99999994 0. ]

[0. 1. ]]

S = [[5.656854 ]

[4.2426405]]

V =

[[ 0.70710677 0.70710677]

[-0.70710677 0.70710677]]

5 eigen结果与手动计算相同

#include 
     
      #include 
      
       #include 
       
          #include 
        
         #include 
         
          #include "opencv2/imgproc/imgproc.hpp"#include 
          
            using namespace std;using namespace cv; //using Eigen::MatrixXf; using namespace Eigen; using namespace Eigen::internal; using namespace Eigen::Architecture; int GetEigenSVD(Mat &Amat, Mat &Umat, Mat &Smat, Mat &Vmat){//-------------------------------svd测试 eigen Matrix2f A; A(0,0)=Amat.at
           
            (0,0); A(0,1)=Amat.at
            
             (0,1); A(1,0)=Amat.at
             
              (1,0); A(1,1)=Amat.at
              
               (1,1); JacobiSVD
               
                 svd(A, ComputeThinU | ComputeThinV ); Matrix2f V = svd.matrixV(), U = svd.matrixU(); Matrix2f S = U.inverse() * A * V.transpose().inverse(); // S = U^-1 * A * VT * -1 //Matrix2f Arestore = U * S * V.transpose(); // printeEigenMat(A ,"/sdcard/220/Ae.txt"); // printeEigenMat(U ,"/sdcard/220/Ue.txt"); // printeEigenMat(S ,"sdcard/220/Se.txt"); // printeEigenMat(V ,"sdcard/220/Ve.txt"); // printeEigenMat(U * S * V.transpose() ,"sdcard/220/U*S*VTe.txt"); Umat.at
                
                 (0,0) = U(0,0); Umat.at
                 
                  (0,1) = U(0,1); Umat.at
                  
                   (1,0) = U(1,0); Umat.at
                   
                    (1,1) = U(1,1); Vmat.at
                    
                     (0,0) = V(0,0); Vmat.at
                     
                      (0,1) = V(0,1); Vmat.at
                      
                       (1,0) = V(1,0); Vmat.at
                       
                        (1,1) = V(1,1); Smat.at
                        
                         (0,0) = S(0,0); Smat.at
                         
                          (0,1) = S(0,1); Smat.at
                          
                           (1,0) = S(1,0); Smat.at
                           
                            (1,1) = S(1,1); // Smat.at
                            
                             (0,0) = S(0,0); // Smat.at
                             
                              (0,1) = S(1,1); //-------------------------------svd测试 eigen return 0;}int main(){ // egin(); // opencv3(); //Eigentest(); //opencv(); //similarityTest(); // double data[2][2] = { { 629.70374, 245.4427 }, // { -334.8119 , 862.1787 } }; double data[2][2] = { { 4, 4 }, { -3, 3} }; int dim = 2; Mat A(dim,dim, CV_64FC1, data); Mat U(dim, dim, CV_64FC1); Mat V(dim, dim, CV_64FC1); Mat S(dim, dim, CV_64FC1); GetEigenSVD(A, U, S, V); Mat Arestore = U * S * V.t(); cout <
                              <