奇异值分解SVD eigen、opencv实现

Theory

先贴上奇异值分解的物理意义：https://www.zhihu.com/question/22237507
其他svd介绍（包括特征分解）：https://www.cnblogs.com/LeftNotEasy/archive/2011/01/19/svd-and-applications.html
https://www.cnblogs.com/pinard/p/6251584.html
再贴一下之前写的笔记：

物理意义参考代码：

#include <opencv2/opencv.hpp>
using namespace cv;

int main(int argc, char** argv)
{
	Mat img = imread("../../cameraman.jpg");
	if (!img.dims)
	{
		return -1;
	}
	int channels = img.channels();
	//分离通道
	std::vector<Mat> src;
	split(img, src);
	for (int i = 0; i < src.size()/*==channels*/; i++)
	{
		src[i].convertTo(src[i], CV_32F);
	}
	//SVD
	std::vector<Mat>W, U, VT,S;
	for (int i = 0; i < channels; i++)
	{
		Mat w, u, vt;
		SVD svd;
		svd.compute(src[i], w, u, vt);
		W.push_back(w);
		U.push_back(u);
		VT.push_back(vt);
	}
	//
	Mat s1 = Mat::zeros(W[0].rows, W[0].rows, CV_32F);
	Mat s2 = Mat::zeros(W[1].rows, W[1].rows, CV_32F);
	Mat s3 = Mat::zeros(W[2].rows, W[2].rows, CV_32F);
	S.push_back(s1); S.push_back(s2); S.push_back(s3);
	for (int j = 0; j < img.rows; j++)
	{
		printf("%d/%d\n", j, img.rows);
		for (int i = 0; i < channels; i++)
		{
			src[i] = Mat::zeros(img.rows, img.cols, CV_32F);
			S[i].at<float>(j, j) = W[i].at<float>(j);
			src[i] = U[i] * S[i] * VT[i];
			src[i].convertTo(src[i], CV_8U);
		}
		merge(src, img);
		imshow("img", img);
		//
		if (j<img.rows*0.05)
			waitKey(1000);
		else
			waitKey(1);
		//
	}
	getchar();
	//
	return 0;
}

eigen实现

方阵SVD分解

#include <iostream>  
#include <Eigen/SVD>  
#include <Eigen/Core>  
#include <Eigen/Dense>    

//using Eigen::MatrixXf;    
using namespace Eigen;    
using namespace Eigen::internal;    
using namespace Eigen::Architecture;    

int main()  
{  
    Matrix3f A;  
    A(0,0)=1,A(0,1)=0,A(0,2)=1;  
    A(1,0)=0,A(1,1)=1,A(1,2)=1;  
    A(2,0)=0,A(2,1)=0,A(2,2)=0;  
    // A(0,0)=1,A(0,1)=1,A(0,2)=1;  
    // A(1,0)=1,A(1,1)=1,A(1,2)=1;  
    // A(2,0)=1,A(2,1)=1,A(2,2)=1; 
    MatrixXf C(3, 1);
    JacobiSVD<Eigen::MatrixXf> svd(A, ComputeThinU | ComputeThinV );  
    Matrix3f V = svd.matrixV(), U = svd.matrixU();  
    C = svd.singularValues();
    Matrix3f  S = U.inverse() * A * V.transpose().inverse(); // S = U^-1 * A * VT * -1  
    std::cout<<"A :\n"<<A<<std::endl;  
    std::cout<<"C :\n"<<C<<std::endl; 
    std::cout<<"U :\n"<<U<<std::endl;  
    std::cout<<"S :\n"<<S<<std::endl;  
    std::cout<<"V :\n"<<V<<std::endl;  
    std::cout<<"U * S * VT :\n"<<U * S * V.transpose()<<std::endl;  
    system("pause");  
    return 0;  
}

运行结果：

非方阵求伪逆

#include <iostream>
#include <Eigen/SVD>
#include <Eigen/Core>

using namespace std;


// 利用Eigen库，采用SVD分解的方法求解矩阵伪逆，默认误差er为0
Eigen::MatrixXd pinv_eigen_based(Eigen::MatrixXd & origin, const float er = 0) {
    // 进行svd分解
    Eigen::JacobiSVD<Eigen::MatrixXd> svd_holder(origin,
                                                 Eigen::ComputeThinU |
                                                 Eigen::ComputeThinV);
    // 构建SVD分解结果
    Eigen::MatrixXd U = svd_holder.matrixU();
    Eigen::MatrixXd V = svd_holder.matrixV();
    Eigen::MatrixXd D = svd_holder.singularValues();
    //cout<<"D :\n"<<D<<endl; 
    // 构建S矩阵
    Eigen::MatrixXd S(V.cols(), U.cols());
    S.setZero();

    for (unsigned int i = 0; i < D.size(); ++i) {

        if (D(i, 0) > er) {
            S(i, i) = 1 / D(i, 0);
        } else {
            S(i, i) = 0;
        }
    }

    // pinv_matrix = V * S * U^T
    return V * S * U.transpose();
}



int main() {
    // 设置矩阵行数、列数
    const int ROW = 3;
    const int COL = 4;

    // 生成大小 ROW * COL 的随机矩阵
    Eigen::MatrixXd A;
    A = Eigen::MatrixXd::Random(ROW, COL);
    //试了下指定元素
    A(0,0)=1,A(0,1)=0,A(0,2)=1,A(0,3)=1;  
    A(1,0)=0,A(1,1)=1,A(1,2)=1,A(1,3)=0;  
    A(2,0)=0,A(2,1)=0,A(2,2)=0,A(2,3)=1;     
    // 打印矩阵A
    cout << "矩阵A为：" << endl;
    cout << A << endl;
    
    // 打印矩阵A的伪逆矩阵
    cout << "矩阵A的伪逆为：" << endl;
    cout << pinv_eigen_based(A) << endl;
}

运行结果：

opencv 实现

#include <opencv2/core/core.hpp>
// #include <opencv2/highgui/highgui.hpp>
// #include <opencv2/imgproc/imgproc.hpp>
#include <iostream>
 
using namespace std;
using namespace cv;


//发现没啥用
void printMat(Mat & m){

	for (int row = 0; row < m.rows; row++){
        for (int col = 0; col < m.cols; col++){
            cout<<m.at<float>(row,col)<<"  ";
        }
		cout<<endl;
	}
}

int main(){

    const int ROW = 3;
    const int COL = 3;
    //Mat mat = Mat::ones(ROW, COL, CV_32F);
    Mat mat = (Mat_<float>(3, 3) << 1,0,1,0,1,1,0,0,1);
    
    cout<<"mat = "<<endl;
    printMat(mat);
    SVD svd(mat);
    Mat UT,V;

    Mat U = svd.u;
    transpose(U,UT);
    cout<<"U = "<<endl;
    cout<<U<<endl;

    transpose(svd.vt,V);
    cout<<"V = "<<endl;
    cout<<V<<endl;

    Mat D = svd.w;
    cout<<"D = "<<endl;
    cout<<D<<endl;
    //printMat(D);
    

    Mat I3 = Mat::ones(3,1,CV_32F);
    Mat D_inv;
    divide(I3,D,D_inv);
    Mat W = Mat::diag(D_inv);
    cout<<"W = "<<endl;
    printMat(W);

    Mat mat_inv = V*W*UT;
    cout<<"mat inverse by SVD :"<<endl;
    printMat(mat_inv);

    Mat mat_inv1 = mat.inv(DECOMP_CHOLESKY);
    cout<<"mat inverse by inv() :"<<endl;
    printMat(mat_inv1);
    return 0;
}

代码没整理，见谅。运行结果：