参考
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 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:
U:
-1.000 0.000 0.000 1.000
S:
5.657 0.000 0.000 4.243
VT
-0.707 -0.707 -0.707 0.707