numpy opencv matlab eigen SVD结果对比

2022-10-24 10:34:48

参考

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

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

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

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

https://*.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://*.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 unitary matrix, 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 手动计算

给定一个大小为的矩阵，虽然这个矩阵是方阵，但却不是对称矩阵，我们来看看它的奇异值分解是怎样的。

由进行对称对角化分解，得到特征值为，，相应地，特征向量为，；由进行对称对角化分解，得到特征值为，，

当 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)】可见特征向量不唯一

特征向量为，。取，则矩阵的奇异值分解为

2 MATLAB 结果与手动计算不同

AB =  [[ 4 4 ],

  [-3   3 ]]

[U,S,V] = svd(AB);

U

S

V'#需要转置

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 np

import math

import cv2

a = np.array([[4,4],[-3,3]])

# a = np.random.rand(2,2) * 10

print(a)

u, d, v = np.linalg.svd(a)

print(u)

print(d)

print(v)#不需要转置

A = [[ 4 4]
[-3 3]]

U =
[[-1. 0.]
[ 0. 1.]]

S=
[5.65685425 4.24264069]

V=
[[-0.70710678 -0.70710678]
[-0.70710678 0.70710678]]

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

import numpy as np

import math

import cv2

a = np.array([[4,4],[-3,3]],dtype=np.float32)

# a = np.random.rand(2,2) * 10

print(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 <iostream>

#include <Eigen/SVD>

#include <Eigen/Dense>

#include <opencv2/core/core.hpp>

#include <opencv2/highgui/highgui.hpp>

#include "opencv2/imgproc/imgproc.hpp"

#include <iostream>

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<double>(0,0);

    A(0,1)=Amat.at<double>(0,1);

    A(1,0)=Amat.at<double>(1,0);

    A(1,1)=Amat.at<double>(1,1);

    JacobiSVD<Eigen::MatrixXf> 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<double>(0,0) = U(0,0);

    Umat.at<double>(0,1) = U(0,1);

    Umat.at<double>(1,0) = U(1,0);

    Umat.at<double>(1,1) = U(1,1);

    Vmat.at<double>(0,0) = V(0,0);

    Vmat.at<double>(0,1) = V(0,1);

    Vmat.at<double>(1,0) = V(1,0);

    Vmat.at<double>(1,1) = V(1,1);

    Smat.at<double>(0,0) = S(0,0);

    Smat.at<double>(0,1) = S(0,1);

    Smat.at<double>(1,0) = S(1,0);

    Smat.at<double>(1,1) = S(1,1);

    // Smat.at<double>(0,0) = S(0,0);

    // Smat.at<double>(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 <<A<<endl;

    cout <<Arestore<< endl;

    cout <<"U " << U<<endl;

    cout <<"S " <<S<<endl;

    cout <<"V " << V.t()<<endl;

    cout <<V<<endl;

    return 0;

}

[4, 4;
-3, 3]

U =

[0.9999999403953552, 0;
0, 0.9999999403953552]
S =

[5.656854629516602, 0;
0, 4.242640972137451]
V =

[0.7071067690849304, 0.7071067690849304;
-0.7071067690849304, 0.7071067690849304]

6 在线计算网站与手动计算不同

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

Input matrix:

 4.000  4.000

-3.000  3.000

Singular Value Decomposition:

-1.000  0.000

 0.000  1.000

5.657 0.000

0.000 4.243

V^T

-0.707 -0.707

-0.707  0.707

码农公寓

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

numpy.linalg.svd 源码

Input matrix:

Singular Value Decomposition:

相关文章